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Abstract 

The cold tandem rolling of metal strip presents a 
significant control challenge because of nonlinearities 
and process complexities. Flatness and edge drop are two 
significant parts of cold rolling strip quality. The flatness 
and edge drop control system is a coupled multivariable 
and nonlinear control system. Analyzing the coupling 
relationship and realizing the decoupling control between 
them are very important parts of achieving the high- 
precision control. In this paper, based on the flatness and 
edge drop control capabilities analysis of all the control 
actuators, the simplified model of flatness and edge drop 
is established, and the coupling relationship between 
flatness and edge drop is analyzed. To control the flatness 
and edge drop precisely in cold rolling mill, the 
decoupling control of flatness and edge drop is very 
important. Inverse decoupling control realizes the 
decoupling control of the original system by connecting 
the inverse system in series before the original system 
and constituting the pseudo linear composite system. It 
has been applied in many industrial processes 
successfully. Combining the neural network which has 
the powerful ability in nonlinear mapping and self- 
learning with inverse system, the neural network inverse 
decoupling control method is proposed to realize the 
decoupling control of flatness and edge drop in this paper. 
Simulation results represent that the flatness and edge 
drop control performance is improved effectively. 

Keywords: Flatness; Edge Drop; Control capability; 
Decoupling Control; Cold Rolling Mill 

1. Introduction 

The cold tandem rolling of metal strip presents a 
significant control challenge because of nonlinearities 
and process complexities. Flatness and edge drop are two 
significant parts of cold rolling strip quality. Nowadays, 
Flatness control precision has satisfied the demand of 
users, but edge drop control level still has no 
breakthrough [1] . The edge drop of strip will directly 
influence the side scrapping of strip and reduce the 
product yield [2] . In recent years, the edge drop control 



has become the concern issue in the rolling process 
control field. The technologies such as K-WRS 
(Kawasaki- Work Roll Shifting) work roll, EDC (Edge 
Drop Control) work roll etc. are proposed one after 
another [3, 4] . However, the control system has not 
achieved good control effect yet. There are a lot of 
problems to resolve for the edge drop control. The edge 
drop is mainly controlled by work roll shifting. Yet, the 
work roll shifting inevitably influences the shape of 
loaded roll gap. Due to the inseparable relations between 
the roll gap and flatness, the flatness must be affected. 
The coupling relationship between the edge drop and 
flatness has been a hot topic [5] . 

The composition control of flatness and edge drop is a 
coupled multivariable and nonlinear control system 
because of the nonlinearities and complexities in cold 
rolling process. So realizing the decoupling control of 
flatness and edge drop is one of the most important parts 
of achieving the high-precision control. At present, the 
edge drop control is compensated by adding a work roll 
bending compensation part during the flatness control, 
the control performance was improved in a certain extent 
[6] . But the decoupling control of flatness and edge drop is 
still not realized in essence. 

Inverse decoupling control realizes the decoupling 
control of the original system by connecting the inverse 
system in series before the original system and 
constituting the pseudo linear composite system. It has 
been applied in many industrial processes successfully. [7] 
Combining the neural network which has the powerful 
ability in nonlinear mapping and self-learning [8, 9] with 
inverse system, the neural network inverse decoupling 
control will gain the better decoupling control effect. 

In this paper, the flatness and edge drop control abilities 
of all the control actuators are analyzed. The coupling 
relationship of flatness and edge drop is analyzed. The 
neural network inverse decoupling control method is 
proposed to realize the decoupling control of flatness and 
edge drop. Simulation results represent that the flatness 
and edge drop control performance is improved 
effectively. 
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The remainder of the paper is organized as follows: 
Section (2) focuses on the establishment of the analytic 
finite element model and the choice of the rolling 
simulation conditions. The elastic deformation of the roll 
and plastic deformation of the strip are considered as a 
whole to analyze the strip rolling process thoroughly. 
Section (3) emphasizes on the work roll bending, shifting 
and intermediate roll bending and shifting control 
capability analysis of the flatness and edge drop. Section 
(4) develops the coupling analysis between flatness and 
edge drop based on the results of section (3). Section (5) 
focuses on the decoupling control of flatness and edge 
drop based on the inverse system and neural network. 
And section (6) is conclusion of this paper. 

2. The establishment of the analytic model 

To deeply analyze the control characteristic and the 
coupling relation of the flatness and edge drop, the 
analytic finite element model is established in which 
elastic deformation of the roll and plastic deformation of 
the strip are considered as a whole. The simulation is for 
five stands 6-high tandem cold rolling mill. The roll finite 
element deformation model is shown as Figure 1. 




Figure 1. Roll deformation model 

The diameter of the backup roll is 1276mm. The diameter 
of the intermediate roll is 467mm. The diameter of the 
work roll is 426mm. The rolling velocity is 1057m/min. 
And the entry tension is 60Mpa, exit tension is 150 Mpa, 
friction coefficient is 0.05, reduction ratio is 0.05, strip 
width is 1100mm, strip thickness is 2.2mm. The other 
simulation factors about the control actuators are shown 
as Table 1 as follows. 



Table 1 The rolling simulation conditions 



Control actuators 


Min 


Avg 


Max 


Bending force of work 
roll (KN/side) 


-180 


0 


360 


Bending force of 
intermediate roll (KN/ 




0 


500 


side) 








Shift of work roll (mm) 


-100 


0 


25 


Shift of intermediate roll 
(mm) 


-100 


0 


50 



aims at the edge of strip; the shift is positive when the 
work roll taper corner enters into the edge of strip and 
negative when the work roll taper comer stretch out the 
edge of strip. The definition of intermediate roll shift is 
as follows: the shift of intermediate roll Sj equals zero 

when the intermediate roll corner aims at the edge of strip; 
the shift is positive when the intermediate roll corner 
enters into the edge of strip and negative when the 
intermediate roll corner stretch out the edge of strip. 
Figure 2 shows the definition of the work roll shift 
S w and the intermediate roll shift S { . 











Figure 2. The definition of the roll shift 

The control of flatness is actually the control of the 
loaded roll gap in the rolling process. The strip cross 
section is mainly symmetrical. The loaded roll gap can be 
shown as quartic polynomial as Equation 1. 

/(x) = a 0 +a 2 x 2 +ti 4 i 4 ,iG [-l,+l] (1) 

Where X is the corresponding width in view of 
considering the roll gap center as the original 
point, a 0 ,a 2 , a 4 are coefficients. 

The roll gap can be divided into quadratic part f 2 (x) , 
quartic part f 4 (x) and constant / 0 (x) . The constant 
/ 0 (x) is expressed as Equation 2. 

fo( X ) = f( ±l ) = a 0 +a 2 +a 4 (2) 

The quadratic part / 2 (x) is expressed as Equation 3. 

f 2 {x) = C w2 (\-x 2 ) (3) 

The curve of the quadratic part f 2 (x) is shown as 

Figure 3. It is related to the quadratic crown and the 
quadratic flatness of the strip. 

The quartic part f 4 (x) is expressed as Equation 4 

f,{x) = f (x)- f 0 (x)- f 2 {x) 

= a 4 ft 4 - x 2 j 



The definition of work roll shift is as follows: the shift of 
work roll S w equals zero when the work roll taper corner 
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The curve of the quartic part f 4 (x) is shown as Figure 

4. It is related to the quartic crown and the quartic 
flatness of the strip. 




Figure3. The quadratic crown of the strip 



crown and edge drop control capabilities of control 
actuators in this paper. 

3. The control capability analysis of the 
flatness and edge drop 

The main control actuators are the work roll bending, 
intermediate roll bending, work roll shifting and 
intermediate roll shifting on 5 stands 6-high cold tandem 
rolling mill. So the quadratic crown, quartic crown and 
edge drop control capabilities of these four control 
actuators are analyzed as follows. 

Figure5 shows the quadratic crown, quartic crown and 
edge drop control capabilities of the work roll bending 
when the intermediate roll shifting is 0mm, the 
intermediate roll bending force is OKN, and the work roll 
shifting is -100mm, 0mm, or 25mm. 




Figure4. The quartic crown of the strip 



The quadratic crown C W2 is defined as subtracting the 

edge value (x = ±l) from the center value (x = 0) of 
the loaded roll gap. It can be obtained as Equation 5. 



C W2 = / (0) — / (±!) = — («2 + ) ( 5 ) 

The quartic crown C W4 is defined as subtracting the 
edge value (x = ±l) from the extreme value 
^x = ±V2 / 2 j of the loaded roll gap. It can be obtained as 
Equation 6. 



C W4 — f 4 



S 



-f*( °)=- 



( 6 ) 



So the roll gap can be expressed as Equation 7 as follows. 

/ M = f 0 ( x ) + f 2 (*) + h M 



— @q+(4-C W 4 4 C w4 x 



(7) 



The quadratic component of the loaded roll gap 
corresponds to the quadratic waves of the strip. The 
quartic component corresponds to the quartic waves. 
Therefore, the control of the quadratic waves is mainly 
for the control of the quadratic crown, and the control of 
the quartic waves is mainly for the control of the quartic 
crown. For this reason, the flatness can be approximately 
expressed by the quadratic crown and the quartic crown. 
So the research is mainly on the quadratic crown, quartic 




(a) the quadratic crown control capabilities 




(b) the quartic crown control capabilities 




(c) edge drop control capabilities 



Figure5. The control capabilities of the work roll bending 
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Figure (a) shows the quadratic crown curve changing 
with the work roll bending. Figure (b) shows the quartic 
crown curve changing with the work roll bending. Figure 
(c) shows the edge drop curve changing with the work 
roll bending. The results represent that the quadratic 
crown, the quartic crown and edge drop all decrease with 
the increasing work roll bending force and the changes 
are nearly linear. But the. quadratic crown control 
capability of work roll bending is most effective. The 
control scope of the quadratic crown is about 220um 
under the work roll bending control. The control scope of 
the edge drop is about 50um. And the quartic crown 
control ability of the work roll bending is about 20um. 
So the work roll bending is mainly used to control the 
quadratic crown. It is also a way to control edge drop. 




(a) the quadratic crown control capabilities 




(b) the quartic crown control capabilities 




(c) edge drop control capabilities 



Figure6. The control capabilities of the intermediate roll bending 



Figure6 shows the quadratic crown, quartic crown and 
edge drop control capabilities of the intermediate roll 
bending when the intermediate roll shifting is 0mm, work 
roll bending force is OKN, and work roll shifting is - 
100mm, 0mm, or 25mm. 

Figure (a) shows the quadratic crown curve changing 
with the intermediate roll bending. Figure (b) shows the 
quartic crown curve changing with the intermediate roll 
bending. Figure (c) shows the edge drop curve changing 
with the intermediate roll bending. The results represent 
that the quadratic crown and edge drop all decrease with 
the increasing intermediate roll bending force and the 
changes are nearly linear. But the quartic crown increases 
with the increasing intermediate roll bending force and 
the change is nearly linear. The. quadratic crown control 
capability of intermediate roll bending is about 50um 
under the intermediate roll bending control. And the 
quartic crown control orientation of the intermediate roll 
bending is same as the work roll bending. The control 
scope of the edge drop is about lOum. So the 
intermediate roll bending is seldom used for the control 
of the edge drop. And the quartic crown control ability 
of the intermediate roll bending is about 20um. And the 
quartic crown control orientation of the intermediate roll 
bending is opposite to the work roll bending. So it is 
important to consider the control strategy of the work roll 
bending and the intermediate roll bending to control the 
quadratic and quartic crown. 

Figure7 shows the quadratic crown, quartic crown and 
edge drop control capabilities of the work roll shifting 
when the intermediate roll shifting is 0mm, intermediate 
roll bending force is OKN, and work roll bending is - 
180KN,0KN, or 360KN. 




(a) the quadratic crown control capabilities 
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(b) the quartic crown control capabilities 




(c) edge drop control capabilities 

Figure7. The control capabilities of the work roll shifting 

Figure (a) shows the quadratic crown curve changing 
with the work roll shifting. Figure (b) shows the quartic 
crown curve changing with the work roll shifting. Figure 
(c) shows the edge drop curve changing with the work 
roll shifting. The results represent that the quadratic 
crown, the quartic crown and edge drop all change little 
with the increasing work roll shifting when work roll 
shifting value is smaller than 0mm, but the quadratic 
crown, the quartic crown and edge drop all decrease with 
the increasing work roll shifting when work roll shifting 
value is greater than 0mm. So the edge drop is controlled 
effectively when the work roll shifting is positive. The 
edge drop can decrease 60um when the work roll shifting 
is from 0mm to 25mm. So the work roll shifting is the 
most effective way to control edge drop. At the same 
time , the quadratic crown decreases about 35um and the 
quartic crown decreases about 20um. So the coupling 
relation exits between the edge drop and flatness control. 
Figure8 shows the quadratic crown, quartic crown and 
edge drop control capabilities of the intermediate roll 
shifting when the work roll shifting is 0mm, intermediate 
roll bending force is 0KN, and work roll bending is - 
180KN,0KN, or 360KN. 




-m bending force -180KN 
-"WR bending force OKN 
-"WR bending force 360KN 



-10 L 
-100 



-50 0 

intermediate roll shift(mm) 



(b) the quartic crown control capabilities 




(c) edge drop control capabilities 



Figure8. The control capabilities of the intermediate roll shifting 

Figure (a) shows the quadratic crown curve changing 
with the intermediate roll shifting. Figure (b) shows the 
quartic crown curve changing with the intermediate roll 
shifting. Figure (c) shows the edge drop curve changing 
with the intermediate roll shifting. The results represent 
that the quadratic crown and edge drop both decrease 
slowly with the increasing intermediate roll shifting. But 
the quartic crown increases slowly with the increasing 
intermediate roll shifting. So in the cold rolling control 
system, the intermediate roll shifting is usually set in the 
preset control system and seldom changed in the 
feedback control. 

In the real rolling process, four control actuators are used 
according to the different control strategy to meet 
different demands for kinds of cold rolling mills. The 
main purpose is to develop the control ability of each 
control actuator to the most effective extent. The analysis 
of the control capability of control actuators is very 
important for the control strategy. In view of above 
analysis, the quadratic crown control capability sequence 
of the actuators is: the capability of the work roll 
bending > the capability of the work roll shifting > the 
capability of the intermediate roll bending > the 
capability of the intermediate roll shifting; the quartic 
crown control capability sequence of the actuators is: the 
capability of the work roll shifting > the capability of the 
work roll bending > the capability of the intermediate roll 
bending > the capability of the intermediate roll shifting; 
the edge drop control capability sequence of the actuators 
is: the capability of the work roll shifting > the capability 
of the work roll bending > the capability of the 
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intermediate roll bending > the capability of the 
intermediate roll shifting. 

4. The coupling analysis between flatness 
and edge drop 

The flatness is approximately expressed by the quadratic 
crown and the quartic crown, the quadratic crown is more 
larger than the quartic crown, thus, the coupling model of 
the flatness and edge drop can be simplified as the 
coupling model of the quadratic crown and edge drop. 
From the control capability of the actuators, it can be 
seen that the strongest control capability of the edge drop 
is work roll shifting, secondly is work roll bending. The 
control of the intermediate roll shifting and bending can 
be ignored. The strongest control capability of the 
quadratic crown is work roll bending, secondly is work 
roll shifting. The control of the intermediate roll shifting 
and bending can be ignored. Hence, the coupling model 
of the flatness and edge drop can be simplified as the 
coupling model of the quadratic crown and edge drop 
controlled by the work roll bending and work roll shifting. 
In the actual rolling process, the quadratic crown is 
mainly controlled by work roll bending. Edge drop is 
mainly controlled by work roll shifting. But edge drop is 
influenced during controlling the quadratic crown by 
work roll bending and coupling relation is existed. 
Equally, the quadratic crown is influenced during 
controlling edge drop by work roll shifting. 

The flatness and edge drop control capability of each 
control actuator can be defined as induced variation of 
the quadratic crown, the quartic crown and edge drop per 
change of the control actuator. The flatness and edge 
drop control capabilities of the work roll bending and the 
work roll shifting are as table 2 shown as follows. 

Table2 The flatness and edge drop control capability of the actuators 



Control 

actuators 


Work roll bending 


Work roll shifting 


the quadratic 
crown 


-0.40736 


-0.24779 


Edge drop 


-0.09094 


-0.49904 



In this paper, the coupling coefficient k e2 is defined as 

the influence for the quadratic crown during controlling 
edge drop by work roll shifting as Equation 8. 



k 



el 



K c s 

^W2^W 

K cs 



( 8 ) 



Where K . ,, is the quadratic crown control 

capability of the work roll shifting; K c is the edge 
drop control capability of the work roll shifting. 

The coupling coefficient k 2e is defined as the influence 

for the edge drop during controlling the quadratic crown 
by work roll bending as Equation 9. 



K cf 

U C eMy 

K le ~ ' 



K, 



Cw?F w 



(9) 



Where K c is the edge drop control capability of the 

work roll bending; is the quadratic crown 

control capability of the work roll bending. 

As a result, k e2 = 0.496533 and k 2e = 0.0223242 , 

the coupling relation is strong during controlling the 
quadratic crown by work roll bending and controlling the 
edge drop by work roll shifting. Especially the influence 
for the quadratic crown during controlling edge drop by 
work roll shifting is very large. The control quality 
characteristic is decreased, so the flatness and edge drop 
need to be decoupled in the control process. 

5. The Decoupling Control of Flatness and 
Edge Drop 

Inverse system can achieve the representation from the 
output to the input of the original system. Inverse 
decoupling control realizes the decoupling control of the 
original system by connecting the inverse system in 
series before the original system and constituting the 
pseudo linear composite system. 

If the original system has no precision model, the analytic 
inverse system can not be established. The tandem rolling 
of cold metal strip presents a significant control challenge 
because of nonlinearities and process complexities. Due 
to the powerful ability in nonlinear mapping and self- 
learning of the neural network, the flatness and edge drop 
of the strip control model can be established by the neural 
network. So the neural network inverse decoupling 
control combining the neural network nonlinear mapping 
capability with inverse system decoupling control will 
gain the better decoupling control effect [10, 11] . Neural 
network inverse decoupling control method achieves the 
dynamic function of the inverse system by adding several 
delay factors or integrators to the static neural network. 
The static neural network approximates the static 
nonlinear function and the delay factor or integrators 
reflect the dynamic characters of the system. The whole 
system can be changed to the decoupled normalized 
pseudo linear composite system. It has the linear transmit 
relation by connecting the neural network inverse system 
in series before the original system. 

Flatness and edge drop control actuators model can be 
simplified to second order dynamic system because the 
work roll bending and work roll shift are controlled by 
the hydraulic cylinder and servo valve [12] . The inputs of 
the inverse system are the quadratic crown and edge drop 
which are the outputs of the original system and first 
order, second order derivative of quadratic crown and 
edge drop. The outputs of the inverse system are work 
roll bending and work roll shift which are the inputs of 
the original system. The pseudo linear composite system 
is established by connecting the neural network inverse 
system and original system in series. The system is 
shown as Figure 9. The inputs of this system are the first 
order, second order derivatives of the edge drop and the 
quadratic crown which are connected to static neural 
network inverse system by four integrators. The outputs 
of the edge drop and the quadratic crown can be obtained 
by neural network inverse system and original system. So 
the flatness and edge drop has been decoupled to two 
independent systems. 



6 



ACSE Journal, Volume 14, Issue 1, ISSN 1687-4811, ICGST LLC, Delaware, USA, August 2014 




ri=C e ( 



yi=C e 



Decoupled as 

>=> 



'1*2 — C w 



y2=c w 



Figure 9. The pseudo linear composite system 

The pseudo linear composite system is equivalent to two 
linear integral subsystems which are two 2-order 
subsystem of quadratic crown and edge drop. Then a 
close-loop controller can be designed by using the linear 
system theorem for the pseudo linear composite system. 
The whole control principle of the flatness and edge drop 
neural network inverse decoupling system is shown as 
Figure 10. The control system simulation is accomplished 
by the MATLAB/SIMULINK. 




Figure 10. The decoupling control of the flatness and edge drop 

The edge drop and flatness control curve of the original 
system without decoupling control is shown as Figure 1 1 . 
It can be seen from Figure 11 that work roll shift value 
changes from 0mm to 5mm at 0.5 second during 
controlling edge drop by work roll shift in the original 
system. The edge drop changes from 19um to 5um, the 
quadratic crown also changes from 35um to 29um. 
During controlling the quadratic crown by work roll 
bending, work roll bending value changes from 0KN to 
80KN at 1 second. It can be seen that the quadratic crown 
changes from 29um to -2um, the edge drop also changes 
from 5um to -2um. Simulation results show that the 
strong coupling relation exists between the edge drop and 
the quadratic crown in the original system. So if the 
decoupling control is not considered, the control effect is 
not ideal and the strip quality is not satisfactory. 




(a) the edge drop output 




(b) the quadratic crown output 

Figure 1 1 . The output of the original system 

The edge drop and flatness control curve of the inverse 
decoupling system is as Figure 12 shown. 




(a) edge drop output 




(b) the quadratic crown output 



Figure 12. The output of the decoupling system 
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It can be seen from Figure 12 that the edge drop changes 
from 20um to Oum at 0.5 second during controlling edge 
drop by work roll shift in the pseudo linear composite 
system, but the quadratic crown nearly has no change. 
And during controlling the quadratic crown by work roll 
bending, the quadratic crown changes from 35um to Oum 
at 1 second, the edge drop changes little, less than 3um. 
Simulation results show that the neural network inverse 
system decoupling system has a good dynamic 
decoupling effect between the edge drop and the 
quadratic crown. 

6. Conclusion 

The flatness and edge drop control capabilities of the 
actuators for the cold tandem rolling mill are analyzed in 
this paper. The simplified model of the flatness and edge 
drop is established and the coupling relation between 
flatness and edge drop is quantitatively analyzed. From 
the analysis results, the edge drop is influenced when 
flatness is controlled by work roll bending and flatness is 
influenced when edge drop is controlled by work roll 
shifting. The strong coupling relation exists between the 
flatness and edge drop. The neural network inverse 
decoupling system of the flatness and edge drop is 
established. The decoupling control is completed by 
combining the inverse system with the original system 
into pseudo linear composite system. Simulation results 
show that the neural network inverse decoupling system 
has a good dynamic decoupling effect between edge drop 
and flatness control. 
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Abstract 

In this context of changing and challenging market 
requirements, Gas Insulated Substation(GIS) has 
found a broad range of applications in power systems 
for more than two decades because of its high 
reliability, easy maintenance and small ground space 
requirement etc. SF 6 has been of considerable 
technological interest as an insulation medium in GIS 
because of its superior insulating properties, high 
dielectric strength at relatively low pressure and its 
thermal and chemical stability. SF 6 is generally found 
to be very sensitive to field perturbations such as 
those caused by conductor surface imperfections and 
by conducting particle contaminants. The purpose of 
this paper is to develop techniques, which will 
formulate the basic equations that will govern the 
movement of metallic particles like aluminum, copper 
in a coated as well as uncoated busduct. 

Keywords: Gas Insulated Substation, Busduct, 

Particle, Epoxy 

1. Introduction 

In this context of changing and challenging market 
requirements, Gas Insulated Substation has found a 
broad range of applications in power systems for more 
than two decades because of its high reliability, easy 
maintenance and small ground space requirement etc. 
SF6 has been of considerable technological interest as 
an insulation medium in GIS because of its superior 
insulating properties, high dielectric strength at 
relatively low pressure and its thermal and chemical 
stability. In uniform field the dielectric strength of 
pure SF6 is approximately three times that of nitrogen 
in same conditions. Although, GIS has been in 
operation for several years, some of the problems 
need full attention. These problems include generation 
of over voltages during switching operations like 
enclosure faults and particle contamination. Sulphur 



Hexafluoride (SF6) is generally found to be very 
sensitive to field perturbations such as those caused 
by conductor surface imperfections and by conducting 
particle contaminants. A study of CIGRE group 
suggests that 20% of failures in Gas Insulated 
Substations (GIS) are due to the existence of various 
metallic contaminations in the form of loose 
particles. [1] However, in non uniform fields it 
depends on many factors such as electrode geometry, 
voltage wave shape, polarity, gas pressure etc. 
although in GIS, a highly non- uniform field would be 
generally avoided, such divergent fields can 
occasionally exist due to electrode surface roughness 
or dust or conducting particles between electrodes. As 
a result of this, breakdown could occur due to local 
field enhancement. 

Several authors conducted experiments on insulating 
particles. Insulating particles are found to have little 
effect on the dielectric behaviour of the gases [1-8]. 
However the presence of atmospheric dust containing 
conducting particles, especially on the cathode, 
reduces the breakdown voltage. 

Conducting particles placed in a uniform ac field lift- 
off at a certain voltage. As the voltage is raised, the 
particles assume a bouncing state reaching a height 
determined by the applied voltage. With a further 
increase in voltage, the bounce height and the corona 
current increase until break down occurs [9]. The lift 
off voltage is independent of the pressure of gas. After 
the onset of bouncing, the offset voltage is 
approximately 30% lower than the lift-off voltage. 
Several methods of conducting particle control and 
de-activation are given. [2]Some of them are: 

a. Electrostatic trapping 

b. Use of adhesive coatings to immobilize the particles 

c. Discharging of conducting particles through 
radiation, and 

d. Coating conducting particles with insulating 
films[3] 
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The work reported in this paper deals with the 
movement of metallic particle in 3 -phase common 
enclosure Gas Insulated Busduct with uncoated and 
coated enclosure. In order to determine the axial 
movement in an uncoated enclosure system, the 
Monte-Carlo technique has been adopted in 
conjunction with motion equation. The present paper 
deals with the computer simulation of particle 
movement in 3 -phase common enclosure GIB with 
coated as well as uncoated system. The specific work 
reported deals with the charge acquired by the particle 
due to macroscopic field at the tip of the particle, the 
force exerted by the field on the particle, drag due to 
viscosity of the gas and random behaviour during the 
movement. Wire like particles of aluminium and 
copper of a fixed geometry in a 3 -phase busduct have 
been considered. The particle movement has been 
calculated in uncoated and coated busducts. The 
movement patterns at various power frequency 
voltages have been obtained. Monte-Carlo technique 
has been adopted for determining the axial movement 
of the particle. It has been assumed that at every time 
step the particle comes out from its location at an 
angle less than 1°. 



field having overcome the forces due to its own 
weight and air drag. The simulation considers several 
parameters e.g. the macroscopic field at the surface of 
the particle, its weight, Reynold’s number, coefficient 
of restitution on its impact to both enclosures and 
viscosity of the gas. During return flight, a new charge 
on the particle is assigned based on the instantaneous 
electric field. 

3. Theoretical Study 

Many authors [1-4] have suggested solutions for the 
motion of a sphere or a wire like metallic particle in 
an isolated busduct system. 

The gravitational force F g acting on a particle of 
mass ‘m’, with ‘g’ the accelerastion due to gravity, is 
given by 

Fg = mg .. (1) 

The charge acquired by a vertical wire particle in 
contact with a bare enclosure can be expressed as [5] : 

It 6 0 1 2 E (t 0 ) •• ( 2 ) 

V “t f 21 A 



2. Modeling Technique 

The Figures 1 & 2 show a typical horizontal three 
phase busduct with and without dielectric coating 
respectively. The enclosures are filled with SF6 gas at 
high pressure. A particle is assumed to be at rest at the 
enclosure surface, just beneath the busbar A, until a 
voltage sufficient enough to lift the particle and move 
in the field is applied. 




A, B, C are the conductors 
Figure 1 . A typical 3-phase common enclosure GIB 




Figure 2. Typical 3 -phase common enclosure GIB with dielectric 
coating 

After acquiring an appropriate charge in the field, the 
particle lifts and begins to move in the direction of 



where Q net is the charge on the particle until the next 
impact with the enclosure, 1 is the particle length, r is 
the particle radius, E(t 0 ) is the ambient electrical field 
at t = t 0 . 

The charge carried by the particle between two 
impacts has been considered constant in the 
simulations. This is a reasonable assumption when the 
applied voltage is low or when the gas pressure is 
high. The electrostatic force 

F e = KQ net E .. (3) 

where K is the correction factor less than unity, 
Q net is the particle charge and E is ambient 
electric field. 



The expression for drag force is given by: 



F d = (t)nr(6nK d (y) + 2.6S6[np g ly{t)}°' S ) ..(4) 



where |LXis the viscosity of the fluid (SF 6 : 15. 5x 10' 6 
kg/m.s at 20°C), r is the particle radius, p g is the gas 
density, 1 is the particle length, K d (y) is a drag 



coefficient. 

The motion equation is given by 

md 2 y 



dt 2 



= F e -mg-F d 



( 5 ) 



Where, y is the direction of motion and F d is the drag 
force. The direction of the drag force is always 
opposed to the direction of motion. For laminar flow 
the drag force component around the hemispherical 
ends of the particle is due to shock and skin friction. 
Very limited publication is available for the 
movement of particle in a three-phase busduct, 
however the equation of motion is considered to be 
same as that of an isolated phase busduct. 
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Monte -Carlo technique 

Solving the above equation gives the radial movement 
of moving metallic particle. But the particle 
geometric configuration of the tips is generally 
uneven and the movement of the particle is not 
unidirectional. So, the particle also moves in axial 
direction because of particle surface roughness at the 
tips and randomness of this axial movement can be 
reasonably simulated by using Monte-Carlo 
technique. For computing axial movement of metallic 
particle, it is assumed that the particle comes out from 
its location at an angle less than 0. At every step of 
movement, A new rectangular random number 
between 0 and 1 is generated at every step of 
movement. Accordingly 0 is modified by multiplying 
it with random number and this fixes the position of 
particle at every time step for determining the new 
location of moving metallic particle in terms of axial 
and radial movements. In the next step again new 
position is computed using particle motion equation 
with newly generated random angles as described 
above. However it is be noted that the usage of Monte- 
Carlo technique yields better axial movements while 
the radial movements remain the same. [4, 10] 

4. Simulation of electric field in 3-phase 
busduct with and without coating 

The charge acquired by the horizontal wire particles 
in contact with a bare enclosure and coated enclosure 
can be expressed as given by Srivastava and Anis [2]. 
The electric field in 3 -phase common enclosure GIB 
electrode system at the position of the particle can be 
written as 

E(t) = E 1 (t) + E 2 (t) + E 3 (t) (6) 

Where, E(t) is the resultant field in vertical direction 
due to the field of three conductors on the surface of 
the particle at the enclosure. 

Ei(t), E 2 (t) and E 3 (t) are the components of the 
electrical field in vertical direction. The gravitational 
force and drag forces are considered as described by 
several authors. 

5. Simulation of Particle Motion 

Computer simulations of the motion of metallic wire 
particles were carried out on GIB of 64mm inner 
diameter for each enclosure and 500mm outer 
diameter with various voltages 245KV,300 KV, 
400KV applied to inner conductors with 120° phase 
difference. 

A conducting particle motion, in an external electric 
field will be subjected to a collective influence of 
several forces. The forces may be divided into 
Electrostatic Force (Fe) 

Gravitational Force (mg) 

Drag Force (Fd) 

Software was developed in C language considering 
the above equations and was used for all simulation 
studies. 



6. Results and Discussions 

Figure 3 shows the movement of aluminium particle 
in radial direction for an applied voltage of 400KV 
rms. During its movement it makes several impacts 
with the enclosure. The highest displacement in radial 
direction during its upward journey is simulated to be 
65mm. As the applied voltage increases the maximum 
radial movement also increases as given in Table 1. 
However, it is noticed that even up to a voltage of 
lOOkv, the particle could not bridge the gap. Further 
calculations may reveal the limiting voltage to enable 
the particle to reach the high voltage conductor. A 
graphical representation of radial movement in 
relation to axial movement is given by Monte-carlo 
technique as shown in Figure 4. The movement of 
copper particle was determined for 400kv with similar 
parameters as above and found to have a maximum 
movement of 18mm in radial direction as shown in 
Figure 5. The movements are also calculated for other 
voltages. 

The movement of copper particles is also given in 
Table 1. It is noticed that the movement of copper 
particle is far less than aluminium particle of identical 
size. This is expected due to higher density of copper 
particle. The axial and radial movements of 
aluminium and copper particles are calculated using 
Monte-Carlo technique for three voltages i.e. 245KV, 
300KV and 400KV with a solid angle of 1°. It is 
significant to note that for all the cases considered, 
there is no change in maximum radial movement, 
even when Monte-Carlo method is applied. A 
relatively high value of axial movement is achieved 
with the random angle of 1°. As expected the axial 
movement of copper particle is lower than the 
aluminium. The movement of aluminium and copper 
particles for a coated enclosure surface is shown in 
Figures 6 and 7 respectively. It can be seen that the 
movement has been far less than that of an uncoated 
system. 



Table 1 : Axial and Radial movement of aluminium and copper 
particles(Simulation time : 1 sec) 



Voltage 


Type 


Max. 

Radial 

Movement 


Monte-Carlo 

(0=i°) 






(mm) 


Radial 

(mm) 


Axial 

(mm) 


245 KV 


A1 


28 


28 


320 


Cu 


13 


13 


280 


300KV 


A1 


32 


32 


340 


Cu 


14 


14 


300 


400KV 


A1 


65 


65 


600 


Cu 


18 


18 


340 
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Table 2: Axial and Radial movement of aluminium and copper 
particles with coated electrodes, thickness 50 pm (Simulation time: 
1 sec) 



Voltage 


Type 


Max. 

Radial 

Movement 

(mm) 


Monte-Carlo 

(0=1°) 






Radial 

(mm) 


Axial 

(mm) 


245 KV 


A1 


0.001 


0.001 


100 


Cu 


0.0001 


0.0001 


90 


300KV 


A1 


0.02 


0.02 


210 


Cu 


0.002 


0.002 


183 


400KV 


A1 


0.4 


0.4 


195 


Cu 


0.014 


0.014 


190 




Figure 3. Particle Movement in 3-phase GIB (400KV/A1 / 0.1mm / 
10mm) 




Axial Movement (mm) 

Figure 4. Particle Movement in 3-phase GIB (400KV/A1 / 0.1mm / 
10mm / 1° Monte-Carlo technique) 




Figure 6 : particle Movement in 3 -phase GIB on coated electrodes 
(400KV/A1 / 0.1mm / 10mm / 50 micro meter) 




Figure 7 : particle Movement in 3 -phase GIB on coated electrodes 
(400KV/Cu / 0.1mm / 10mm / 50 micro meter) 



7. Conclusions 

A model has been formulated to simulate the 
movement of wire like particle in 3 -phase common 
enclosure GIB on bare as well as coated electrodes. 
The results obtained are presented and analyzed. 
Monte-carlo simulation is also adopted to determine 
axial as well as radial movements of particle in the 
busduct. Distance traveled in the radial direction is 
found to be same with or without Monte-carlo 
simulation. The displacement on coated electrodes are 
found to be very less compared to uncoated electrode 
system. 
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Abstract 

In a deregulated power system structure, power 
producers and customers share a common 
transmission network for wheeling power from the 
point of generation to the point of consumption. All 
parties in this open access environment may try to 
produce the energy from the cheaper source for 
greater profit margin, which may lead to overloading 
and congestion of certain corridors of the transmission 
network. This may result in violation of line flow, 
voltage and stability limits and there by undermine the 
system security, utilities therefore need to enhance 
“available transfer capability (ATC)” to ensure the 
system reliability is maintained while serving a wide 
range of bilateral and multilateral transactions. This 
paper focuses on the enhancement ATC using multi 
FACTS devises i.e combination of Thyristor Controlled 
Series Capacitor(TCSC) and Thyristor Controlled Phase 
Angle Regulator (TCPAR). The optimal location of 
FACTS devices were determined based on Sensitivity 
methods. The Reduction of Total System Reactive Power 
Losses Method was used to determine the suitable 
location of TCSC and TCPAR for ATC enhancement. 
The effectiveness of proposed method is demonstrated 
on modified IEEE-9 bus system and 24 bus real time 
system using power world simulator 8.0. 1 

Keywords: Deregulation, Available Transfer Capability , 
TCSC, TCPAR, Multi PACTS devises. 

1. Introduction 

The introduction of deregulation has brought several new 
entities in the electricity market place, while on the other 
hand redefining the scope of activities of many of the 
existing players. With the ongoing expansion and growth 
of the electric utility industry, numerous changes are 
continuously being introduced. 

Improved utilization of the existing power system is 
provided through the application of advanced control 



1 This study has been implemented on Power World 

Simulator at NRI Institute of Technology, Agiripalli. 



technologies. Power electronics based equipment, or 
Flexible AC Transmission Systems (FACTS), provide 
proven technical solutions to address these new operating 
challenges being presented today. Thyristor Controlled 
Series Capacitor (TCSC), Thyristor Controlled Phase 
Angle Regulator (TCPAR), Unified Power Flow 
Controller (UPFC) and Static VAR Compensator (SVC) 
are some of the commonly used FACTS controllers. In 
many deregulated markets, the power transaction between 
buyer and seller is allowed based on calculation of ATC. 
Low ATC signifies that the network is unable to 
accommodate further transaction and hence does not 
promote free competition. FACTS controllers like TCSC, 
TCPAR can help to improve ATC by allowing more 
power transactions. 

Electrical energy plays an important role in today’s 
modern society. It is an efficient energy used widely for 
lighting, communication systems, transportation systems, 
industrial purposes and many other areas. In past, the 
electricity industry is in the control of government and 
therefore monopolistic. However over the past decade, 
the electric industry in many countries has undergone 
significant changes to reform into a free market, which is 
also known as deregulation. 




Figure 1 . Typical structure of a vertically 
integrated utility 



Figure 1. shows the typical structure of a vertically 
integrated utility where links of information flow existed 
only between the generators and the transmission system. 
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Similarly, money (cash) flow was unidirectional, from the 
consumer to the electric utility. One of the first step in the 
restructuring process of the power industry is the 
separation of the transmission activities from the 
electricity generation activities. The subsequent step was 
to introduce competition in generation activities. 

An important point to note is that the restructuring 
process was started with the breaking up of a large 
vertically integrated utility which was characterized by 
the opening up of small municipal monopolies to 
competition. 

A system operator for the whole system was appointed 
and it was entrusted with the responsibility of keeping the 
system in balance, i.e. to ensure that the production and 
imports continuously matched the consumption and 
exports. Quite appropriately, the system operator came to 
be known as the Independent System Operator (ISO). 




Figure 2. Typical structure of a deregulated 
electricity system 



Figure 2 shows the typical structure of a deregulated 
electricity system with links of information and money 
(cash) flow between various players. The possibility of 
having such a complex nature of information flow has 
been one of the driving factors in the process of 
deregulation of the power sector. This has been possible 
due to the rapid developments in the fields of 
communication and information technology during the 
nineties decade. 



II 


S(i) 


iiW 

% 


y m 

1 



Figure 3. TCSC Module 



A TCR includes a pair of anti-parallel thyristors that are 
connected in series with an inductor. In a TCSC, a metal 
oxide varistor (MOV) along with a bypass breaker is 
connected in parallel to the fixed capacitor for 
overvoltage protection. A complete compensation system 
may be made up of several of these modules. 

2.2ThyristorControlledPhaseAngleRegulator 

TCPAR is identical with a phase shifting transformer 
with a thyristor type tap changer and plays an important 
role in increasing load ability of the existing system and 
controlling the congestion in the network. FACTS device 
like TCPAR can be used to regulate the power flow in the 
tie-lines of interconnected power system. 

This is also known as Static Phase Shifter (SPS) and 
phase shift with respect to the bus voltage is achieved by 
adding or subtracting a variable voltage component in 
quadrature with the bus voltage. The quadrature voltage 
is injected in series with the transmission line by a 
boosting transformer. The basic arrangement is shown in 
figure 4. 



2.FACTS CONTROLLERS 
2.1Thyristor Controlled Series Compensator 

Thyristor-controlled series capacitor (TCSC) is a 
capacitive reactance compensator, which consists of a 
series capacitor bank shunted by a thyristor controlled 
reactor in order to provide a smoothly variable series 
capacitive reactance. 

Even though a TCSC in the normal operating range is 
mainly capacitive, but it can also be used in an inductive 
mode. The power flow over a transmission line can be 
increased by controlled series compensation with 
minimum risk of sub synchronous resonance (SSR). 

TCSC is a second generation FACTS controller, which 
controls the impedance of the line in which it is 
connected by varying the firing angle of the thyristors. A 
TCSC module comprises a series fixed capacitor that is 
connected in parallel to a thyristor controlled reactor 
(TCR). 






Vi — U> 



EbcciUng Transformer 



Converter 



rv\s\ 

Boosting Transformer 

v* 

V s 



Figure 4. Shows Static Phase Shifter 



3. MODELING OF FACTS DEVICES 

For enhancing of transfer capability, the static models of 
these controllers are considered. It is assumed that the 
time constants in FACTS devices are very small and 
hence this approximation is justified. 
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3.1 Analysis of Transmission Lines and its 
Power Flows and Loss 

Let the complex voltages at bus i and bus j be denoted as 
Vi L 5i and Vj L 5j respectively. 




Figure 5. Model of a Transmission line. 



The complex power flowing from bus i to bus j can be 
expressed as: 

s;=p,-jaj=v;i,=v;(i R+ i c ) <d 

= KVy-y)(G,+z?,)+y ob c )] 

= y 2 [G, + j(fi 9 + B c )] - vy, (G, + B tj )\ 

= [vfq -vyp tj cos (5. -4) +vy/i si.K4 -4)] 

+B c )-vy ] B tj C o S (4 -fy-vypq sii<4 -4)] 

The Active & Reactive power flow form bus i to bus j is, 

P„ =V?G„ -TO cos^TOj B„ sm(4 ; )] ( 2 ) 

Qy=-^(ByTO + TOyC0s<4,) 

-TO® <44 

Where <5. . = <5 — 8- 

Similarly, Real & Reactive Power flows from bus j to bus 
i can be, 

r =v;G„ -TO “ s <44+'1' / j B,y sin(4 ; .)] <4> 

Q, =-^ 2 <b,+ 4)+0( (4)+Wt- ra <4> < 5 > 

The active and reactive Power loss in the line can be 
calculated as, 



P L = y-G tj - vyp 9 008 ( 4 , ) - vyjBq Sin ( 4 , ) 

+v j g v - yypij co ^j > + v i v J B ‘j sin (^ > 

P L = X 2 Gy + vj G, - 2V i V j G ij cos 4 ; . (6) 

Q l =Q,, +Q ;1 

q l = - y 2 (B,., + B c ) + y.v.B,.. cos (4 . ) - v,. v, c tJ s in (4, ) 
-V 2 (B,., + B c ) + v ; .v, 5, . s/n(4 l7 ) + yv, G ; . . cos(4. ) 

Ql= -y 2 (B l7 +B e ) -y 2 (B l7 +B e )+2yy ^ cos(4 7 ) (7) 

3.2 Power Injection Model Of Thyristor 
Controlled Series Compensator (TCSC) 

Thyristor controlled series compensators (TCSC) are 
connected in series with the lines. The effect of a TCSC 
on the network is as a controllable reactance in the related 
transmission line which compensates the inductive 
reactance of the line results in reducing the transfer 
reactance between the buses. This leads to an increase in 
the maximum power that can be transferred on that line in 
addition to a reduction in the effective reactive power 
losses. The series capacitors also contribute to an 
improvement in the voltage profiles. 




Figure 6 shows a model of a transmission line with a 
TCSC connected between buses i and j. The transmission 
line is represented by its lumped ^-equivalent parameters 
connected between the two buses. During the steady state, 
the TCSC can be considered as a static reactance -jX c . 
This controllable reactance, X c is directly used as the 
control variable to be implemented in the power flow 
equation. 

Let the complex voltages at bus i and bus j be denoted as 
Vi L 5i and Vj L 5j respectively. 

The expressions for real and reactive power flows from 
bus i to bus j can be written as from 

p„ = V*<f v -vyfifj.j cos(4j) + B„ sin(<S2)] (8) 

Qy =-V<By+B4-WG,sin (4 ; )+S # cos(4y) (9) 



Pl = P,i +P,i 



Similarly, the real and reactive power flows from bus j to 
bus i can be expressed as, 



P -V G 

ji j ij 



-vy [G cos(S )- 



B,sin(4 l7 )] 



( 10 ) 
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Qj, ,+BJ+VyiG, sin (S IJ )+B I cos(4,) (n) 

The active and reactive power losses in the line with 
TCSC are 

p\ = />; + p, = g u ( v 2 + v; > - 2 f, v,g ; : . cos s^m 

Ql=Q,1 + Qi, 

= -(V, 2 + V 2 )(Bi,.+ B c ) + 2 vyp u cos 4,. < 13 > 

Line flows 

Without TCSC (- jX c ) 






ij ' J ~ij 2 , 2 J 2 , 2 

^ +-v r u +x u 



WithTCSC(- }X c ) 
G„+jB u = 



1 



* r tJ -Kx tJ -x c ) 



r u + j(x u -x c ) r u -j(x u -x c ) 



where 



G ,=- 



r u 



’■ ij 2 -j(x ij -x c ) 



and 

Hence, AG.. = G. . -G. . (without x c - with :* c ) 






= r . 



,. J +x,.. r u -./(.x 0 .-.Y.) 



AG, 



A 5,,.= 



x^.(x c -2x i .) 



(r 2 i ,+ x 2 1 .,){r 2 i ,-;(x,..-x c ) } 
-x c (^ 2 ~x,/ + x c x..) 
(r 2 ,..+ x 2 ,..){r 2 ,..-;(x ( ..-x c ) 2 } 



(14) 

(15) 



P. =V/AG,-\/V ) |AG, cos(4) + AB, sin(<p (17) 



Q; =-V i 2 AB ij -V l V j [AG ij sm(S ij )-A-B ij cos(S ij )] (18) 

Qj = -V, 2 AB..+ l(y.[AG.. sin($,) + AB y cos(A)] (19) 

These equations are used to structure the TCSC to 
Enhance the Power Transfer Capability. 

3 .3 Power Injection of Thyristor Controlled 
Phase Angle Regulator (TCPAR) 

In a Thyristor controlled phase angle regulator, the phase 
shift is achieved by introducing a variable voltage 
component in perpendicular to the phase voltage of the 
line. The static model of a TCPAR having a complex tap 
ratio of l:aLa and a transmission line between bus i and 
bus j is Shown in Figure5 




Figure 7. Model of TCPAR 



The real and reactive power flows from bus i to bus j can 
be expressed as 

S-j =Py- jQy 

= V* i [a 1 V-a*V j MG lj + B ij \ 

= V> 2 1/.(G 0 . + jB tj ) -ry/ (Gjj+jBjj) 

= V 1 i a\G ij +jB ij )-Vy j a{G ij +jB ij )Z{8 j -8,-a) 

= VV G ij+ jV 2 i a 2 B iJ -V l V J a(G IJ + jB ij )cos(8 J -8 t -a) 
-]Vy j a{G ij +jB ij )s,m{8 j -S t -a) 

= V 2 fl G,. .+ jV 1 i a-B i! -VyjaG^cosiSj -S t -a) 

-Wy ] aB ii cos(8j -S i -a) 

-jV;y.aG,. .sin(A. -S : -a) + jV i V j aB ij sin(S j -8 i -a) 



Hence, the Change in line flows due to Series 
Capacitance, 

Hence, the Change in line flows due to Series 
Capacitance, 

The real Power injection at bus T 

p; =^ 2 Aq-^V.[Aq.cos(4.)+AB i . sin(<p (16) 



= [VW G. j-V'VjaG, ,.cos[-(4 -Sj +a) 
+V i V j aB ij sm[-(S i -S j +a)]] 
+j[V 2 ,a 2 B ij -Vy i aB ij cos[-(S i -S j +a) 
-Vy j aG i y n [-(S i -S j +a)]] 
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= [VW G i -Vy j aG ij cos[S i -S j +a) 
-Vy^i B t j sin [A' ; -Sj + a)\ 

+ )[V 2 i a 2 B ij -VV ! aB ij cos[<3; -8 j + a) 
+Vy j aG i yn[8 i -S j +a)] 



A Pj = -aVV i \G u s in (A,. - A .) - B tj cos (S t - A,)] (29) 

(21) A Q, - ay^B^+aVy^ cos(8 i ~8 j )~ B tj s in ($-*,)] (30) 

A Qj = -aV'VjlG'jCosCS'-S^-BySintf'-Sj)] (31) 



The real and reactive power flows from bus i to bus j are 
Complex power, S tj = R + jQ IJ 

P u =R e[S*,]= Re{V> 2 V;. -aV^} 

p » = irVe, -VypG,! cos(^ -Sj + a) 
-V^aB,.sin(^-5,+a) 



These equations will be used to model the TCPAR in the 
Enhancement of transfer capability 




Figure 8. The injection model of the TCPAR 



••• Qij =-Im[S;.]= Im {V'laW.-aVjWy} 

Q, - -tV 2 a 2 B ij -Vy j aB, cos ($ - S i+ a ) 

+ VjVjd Gy sin (Si-Sj + a)] 

Q, - - l /. 2 a 2 B, + Fl/.aB, cos (S, - S i+ a ) 
-VyjaGg sin (<5)-^ +a)] 

Similarly, the real and reactive power flows from bus j to 
bus i are 

P jl =mV* l [V j -aV l ]Y ij } 

p ji = vl i G ij -Vy j aG ij cos[S i -S j +a) 
+Vy j aB ij sm[8 i -S j +a)\ 

And Q ij =-\m{V\[V j -aV i ]Y ij ) 

Q V ' V J a cos (4 S,+a) 

+ Vy,aG„ sin ($-£.+«)] 

1 J l J 1 J 



(24) 



(25) 



The real and reactive power loss in the line having a 
TCPAR can be expressed as 



P L =P ij +P Jl 



P L = l/.VG, + VjGy - 2 vy s a G :j cos (8 l - S j + a) (26) 

Q L =e, ; +Q ;i 



Q,=-y 2 rrB,.-V;B, ; +2yV ; cos(4 +a) (27) 

This mathematical model makes the Y-bus asymmetrical. 
In order to make the Y-bus symmetrical, the TCPAR can 
be simulated by augmenting the existing line with 
additional power injections at the two buses. The injected 
active and reactive powers at bus i (AP), AQJ and bus j 
(APj, AQjj are given as: 

AP i =-o 2 V i 2 G ij -aVy j [G ij sin(S i -S j )-B ij cos(S i -S j )] (28) 



4.0PTIMAL LOCATION BASED ON 
SENSITIVITY APPROACH FOR TCSC 
AND TCPAR DEVICES 

The static conditions are considered here for the 
placement of FACTS devices in the power system. The 
objectives for device placement may be one of the 
following: 

1 . Reduction in the real power loss of a particular line 

2. Reduction in the total system real power loss 

3. Reduction in the total system reactive power loss 

4. Maximum relief of congestion in the system 



4.1. Reduction of total system VAR 
power loss 

Here, a method based on the sensitivity of the total 
system reactive power loss (Q L ) with respect to the 
control variables of the FACTS devices 
For each of the two devices considered in Section 2.1 and 

2.2, the following describes control parameters 
considered: 

• Net line series reactance (Xy) for a TCSC placed 
between buses i and j, 

• Phase shift (aij) for a TCPAR placed between 
buses i and j. 

The reactive power loss sensitivity factors with respect 
to these control variables may be given as follows: 



1. 



2 . 



Foss sensitivity with respect to control 
parameter Xy of TCSC placed between buses i 
and j, 



oQ l 



a = 

1 Sx. 



(32) 



Foss sensitivity with 
parameter 0/ y of TCPAR 
placed between buses i and j, 

bq l 



respect to control 






80, 



( 33 ) 



19 



ACSE Journal, Volume 14, Issue 1, ISSN 1687-4811, ICGST LLC, Delaware, USA, August 2014 



These factors can be computed for a base case power 
flow solution. Consider a line connected between buses i 
and j and having a net series impedance of X ij? that 
includes the reactance of a TCSC, if present. In that line 
0ij is the net phase shift in the line and includes the effect 
of the TCPAR. The loss sensitivities with respect to Xy 
and By can be computed as: 

f)0 R 2 —X 2 

a u =^=[Vf +Vf -2 vy. cos(4 -£)]— | A (34) 

1 7 1 






(35) 



4.2. Selection of optimal placement of 
FACTS devices 

Using the loss sensitivities as computed in the previous 
section, the criteria for deciding device location might be 
stated as follows: 

1. TCSC must be placed in the line having the most 
positive loss sensitivity index CL - . 

2. TCPAR must be placed in the line having the highest 
absolute value of loss sensitivity index by. 

5. Results and Discussions: 

Case studies are conducted on a 5 number of systems - 
IEEE 9, IEEE 24 rts .bus systems. 

For each system, the enhancements of ATC and voltage 
profile at buses are determined with individual FACTS 
devices placed in turn and multi type FACTS devices. 
The following subsections describe the models and study 
conducted in detail. 

The calculation of ATC is carried out by using Power 
world simulator to compute the power flow of each 
transfer case. The optimal power flow under normal and 
at contingency conditions has been carried out. 

The static models of the FACTS controllers as given in 
Section 2.1, .2.2 are considered. TCSC is represented as a 
static reactance and TCPAR is considered as a 
transformer with a complex tap ratio. The optimal 
locations for placing each of these devices are determined 
by sensitivity analysis. 



5.1. Case Study -1-IEEE 9 Bus system 

5.1.1. System Model 

This system consists of 9 buses, 9 line sections, 3 
generator buses and 3 load buses, see figure 8. 



line reactance and TCPAR is placed in line 3(3-6) with a 
phase shift of 2.9 degrees respectively. 




Figure 9. Single line diagram of IEEE 9 bus 
System 



Line 

No(k) 


From-To 

Bus 


Loss Sensitivities 


TCSC(ajj) 


TCPAR 

(by) 


1 


1-4 


-0.5771 


1.441 


2 


8-2 


- 2.765 


-3.26 


3 


3-6 


-0.721 


1.696 


4 


4-5 


-0.095 


0.588 


5 


9-4 


-0.258 


-0.752 


6 


5-6 


-0.334 


-1.19 


7 


6-7 


-0.0785 


0.4486 


8 


7-8 


-0.58 


-1.513 


9 


8-9 


-0.711 


1.679 



5.1.2. Sensitivity Index 

For this system, from Table 1, three cases have been 
considered 

1. A TCSC is placed in lines 4(4-5) and 7(6-7), operated 
with an inductive reactance of 75% and 20% of the line 
reactance; 

2. A TCPAR is placed in the lines 3(3-6), 9(8-9) operated 
with a phase shift of 2.9, 4.5 degrees and unity tap ratio; 

3. For Multi type FACTS (i.e. combination of TCSC & 
TCPAR), TCSC is placed in line 7(6-7) with 20% of the 



TABLE 1. VAR loss sensitivity index of 
IEEE-9 Bus System 
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5.1.3. ATC Enhancement 



Fr 

om 

To 

Bu 

s 


ATC (MW) 


With 

out 

FAC 

TS 


With 

TCS 

C 


% 

Enh 

ance 

men 

t 


With 

TCP 

AR 


% 

Enha 

ncem 

ent 


With 

TCS 

C- 

TCP 

AR 


% 

Enh 

ance 

ment 


1-4 


170.4 


170.4 


N.E 


170.9 


0.251 


171.4 


0.583 


2-8 


124 


124 


N.E 


124 


N.E 


124 


N.E 


3-5 


112.9 


121.6 


7.103 


113.0 


0.044 


122.9 


8.13 


3-8 


135.6 


135.9 


0.228 


142.6 


4.920 


145.6 


6.865 



Table 2. Available Transfer Capability for an 
IEEE 9-Bus System [without Contingency] 

It is observed from the Table 2, the enhancement of ATC 
is of 8.13% with Multi type FACTS which is more than 
that due to placement of individual FACTS devices. 



0.98 



0.97746 



0.97 



0.96 



0.95 



0.94 



a97585 - 0.97444 °' 97585 

[p at.iH [ ^^ t18 






Without FACTS WithTCSC WithTCPAR WithTCSC-TCPAR 
I Without Contingency I With Contingency 



Figure 8. Comparison of voltages with and 
without FACTS under normal and contingency 
condition at bus 5 



Fro 
m - 
To 
Bus 


ATC (MW) 


With 

out 

FAC 

TS 


Wit 

h 

TCS 

C 


% 

Enh 

ance 

men 

t 


Wit 

h 

TCP 

AR 


% 

Enh 

ance 

men 

t 


Wit 

h 

TCS 

C- 

TCP 


% 

Enh 

ance 

men 

t 


1-4 


169.95 


170.0 


0.058 


171.3 


0.822 


172.9 


1.74 


2-8 


87 


87 


N.E 


87 


N.E 


87 


N.E 


3-5 


166.45 


166.5 


0.036 


166.5 


0.048 


167.4 


0.591 


3-8 


166.45 


166.5 


0.036 


166.5 


0.048 


167.4 


0.591 



Table 3: Available Transfer Capability for an 
IEEE 9-Bus System [with Contingency] 



It is observed from the Table 3, the percentage 
enhancement of ATC is of 1.740% with Multi type 
FACTS which is more than that due to placement of 
individual FACTS devices under contingency condition 
at line 8 (i.e line between bus 7 and bus 8 disconnected). 

5.1.4. Voltage Profiles 

It is also observed from the Fig 8 that with multi-type 
FACTS devices the voltage at bus 5 is increased from 
0.96064p.u to 0.9618 p.u under contingency conditions 
compared to individual FACTS devices 



5.2. Case Study -2-IEEE 24 rts Bus system 
5.2.1. System Model: 

This system consists of 24 buses, 38 line sections, 11 
generator buses and 17 load buses. 




Figure 9. Single line diagram of IEEE 24 rts bus 
System 
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5.2.2. Sensitivity Index 



Line 

No(k) 


From-To 

Bus 

i-j 


Loss sensitivities 


TCSC(aij) 


TCPAR(bij) 


1 


1-2 


-0.01232 


0.2421 


2 


1-3 


-0.0474 


-0.256 


3 


1-5 


-0.2897 


1.227 


4 


2-4 


-0.154 


0.6933 


5 


2-6 


-0.3001 


1.149 


6 


3-9 


-0.0756 


0.5534 


7 


3-24 


-4.898 


-4.33 


8 


4-9 


-0.1096 


-0.7256 


9 


5-10 


-0.06604 


-0.129 


10 


6-10 


-1.0578 


-2.146 


11 


7-8 


-1.1189 


2.218 


12 


8-9 


-0.1122 


-0.738 


13 


8-10 


-0.116 


-0.3062 


14 


9-11 


-1.164 


-2.1447 


15 


9-12 


-0.01225 


-2.44029 


16 


10-11 


-2.922 


-3.1667 


17 


10-12 


-3.176 


-3.4889 


18 


11-13 


-0.877 


-1.638 


19 


11-14 


-3.322 


-3.614 


20 


12-13 


-0.368 


-1.185 


21 


12-23 


-4.986 


-4.608 


22 


13-23 


-4.72 


-4.5339 


23 


14-16 


-13.67 


-7.308 


24 


15-16 


-1.3055 


2.351 


25 


15-21 


-4.47 


-4.2029 


26 


15-21 


-4.47 


-4.2029 


27 


15-24 


-4.575 


4.151 


28 


16-17 


-9.865 


-6.384 


29 


16-19 


-1.409 


2.421 


30 


17-18 


-3.43 


-3.594 


31 


17-22 


-1.76 


-2.825 


32 


18-21 


-0.328 


-1.233 


33 


18-21 


-0.328 


-1.233 


34 


19-20 


-0.214 


-0.571 


35 


19-20 


-0.214 


-0.571 


36 


20-23 


-1.0022 


-1.8325 


37 


20-23 


-1.0022 


-1.8325 


38 


21-22 


-2.218 


-3.201 



Table 4. VAR loss sensitivity index of IEEE-24 
rts Bus System 



For this system, from Table 4, three cases have been 

considered (from Section 3.4.2) 

1. A TCSC is placed in lines 1(1-2) and 15(9-12), 
operated with an inductive reactance of 75% and 
20% of the line reactance; 

2. A TCPAR is placed in the lines 23(14-16), 28(16-17) 
operated with a phase shift of 2.5, 4.5 degrees and 
unity tap ratio; 

3. For Multi type FACTS (i.e. combination of TCSC & 
TCPAR), TCSC is placed in line 15(9-12) with 20% 
of the line reactance and TCPAR is placed in line 
23(14-16) with a phase shift of 2.5 degrees 
respectively. 



5.2.3. ATC Enhancement 
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AR 




1-5 


155.0 


155.8 


0.50 


157.9 


1.81 


158.3 


2.046 


23 - 
















15 


817.2 


821.2 


0.49 


817.2 


N.E 


827.2 


1.208 


13 - 
















21 


917.9 


938.6 


2.20 


927.6 


1.04 


947.6 


3.135 


16 - 
















21 


1118.3 


1123.8 


0.49 


1118.3 


N.E 


1183.5 


5.511 


22 - 
















19 


332.4 


332.4 


N.E 


332.4 


N.E 


337.4 


1.493 



Table 5. Available Transfer Capability for an 
IEEE 24 rts Bus System [without Contingency] 



It is observed from the Table .5, the Enhancement of 
ATC is of 5.511% with Multi type FACTS which is 
more than that due to placement of individual FACTS 
devices. 

It is observed from the Table 6, the percentage 
enhancement of ATC is of 5.66% with Multi type 
FACTS which is more than that due to placement of 
individual FACTS devices under contingency condition 
at line 10 (i.e line between bus 6 and bus 10 
disconnected). 
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0 

B 
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TCP A 
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Enh 
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nt 


With 

TCS 

C- 

TCP 

AR 


% 

Enh 

ance 

men 
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1- 

5 


212.31 


212.31 


N.E 


215.5 


1.48 


224.2 


5.32 


23 

15 


846.5 


846.5 


N.E 


846.5 


N.E 


862.5 


1.86 


13 

21 


1007.7 


1101.0 


8.47 


1037.5 


2.87 


1068.2 


5.66 


16 

21 


1159.7 


1185.0 


2.13 


1159.7 


N.E 


1228.1 


5.56 


22 

19 


322.66 


322.66 


N.E 


322.66 


N.E 


326.9 


1.31 



Table 6. Available Transfer Capability for an 
IEEE 24 rts Bus System [with Contingency at 
line 10 (6-10) 



5.2.4. Voltage Profiles 



1 

0.99 
0.98 
0.97 
0.96 
0.95 
0.94 

I Without Conti ngency I With Contingency 



0.97892 



0.97302 ■ 0 . 974 ! 

■ 096731 



0.9614 




0 : 97 5 3 5 — 

0.9743 

I 





III! 



Without FACTS WithTCSC WithTCPAR WithTCSC-TCPAR 



Figure 10. Comparision of voltages with and 
without FACTS under normal and contingency 
condition at bus 24 



It is observed from the Fig 10 that with multi-type FACTS 
devices the voltage at bus 24 is increased from 
0.97302p.u to 0.98461 p.u under contingency conditions 
compared to individual FACTS devices. 



6. Conclusions 

It is observed that the enhancement of Available Transfer 
Capability comparison between various buses of an 
IEEE-9, 24rts, bus systems. Also observed that the 
enhancement of ATC using Multi type FACTS devices is 
more than individual FACTS devices like TCSC, TCPAR 
by placing optimally with sensitivity approach. 

It is observed that the enhancement of ATC is not only 
for IEEE test systems but also for Indian utility system 
with and without FACTS devices under normal and 
contingency conditions. 

Also observed the voltages with and without FACTS 
devices while enhancing the Available Transfer 
Capability. This indicates that the system is more voltage 
stable while enhancing ATC with multi type FACTS 
devices than with individual FACTS devices under 
contingency conditions. 

Nomenclature 

Py Real power flow from bus i to j 
Qij Reactive power flow from i to j 
P L Real power losses in the line. 

Q l Reactive power losses in the line. 

Pi 1 Real power injection at bus i 

Qi 1 Reactive power injection at bus i 

Pij 1 Real power flow from bus i to j with TCSC 

Qij 1 Reactive power flow from bus i to j with TCSC 

Pl 1 Real power losses in the line with TCSC 

Ql 1 Reactive power losses in the line with TCSC 

Vi L 5i Complex voltage at bus i 

Vj L 8 Complex voltage at bus j 
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Abstract 

In this paper, a reduced order observer is designed and 
implemented for state estimation of an internal 
combustion engine (ICE) model. In most engineering 
systems, certain states are directly measurable. Therefore, 
it is not cost effective and increases computational 
complexity to use a full order observer which will 
demand faster, more powerful processors. This is 
particularly the case in powertrain systems such as 
internal combustion engine. Here, we validate the use of 
a reduced order observer as opposed to a full order 
observer as a suitable solution to this problem. This is 
done through comparison, simulation and analysis of the 
results. 

Keywords: Reduced order observer, Internal combustion 
engine, Linear systems modelling, State estimation. 

1. Introduction 

The introduction of a reduced dimension observer for the 
case where some of the state information is readily 
available was first introduced by Luenberger in [1]. The 
impact of the reduced order observer in state estimation 
and control has been greatly valuable ever since. 
Recently, various control strategies and schemes have 
been designed and implemented using a reduced order 
observer for state and/or parameter estimation [2-7]. 
Such systems include microalternators, induction motors, 
servo motors, reluctance motors, automatic transmissions 
and permanent magnet synchronous machines [8-12]. 
Previous work involving the design and implementation 
of a reduced observer can be found in this journal’s past 
publications [13-14]. 

In all of the aforementioned systems, the reduced order 
observer was used on linear systems with directly 
measurable states to estimate states that are 
immeasurable otherwise. The advantage of reducing the 
observer’s dimension is obvious; a lower order observer 
is analytically less complex. In practice, this means 
fewer sensors, which reduce the cost of the overall 
system. Reduced order ICE models can be found 
throughout literature [15-17]. Similarly to the work 



presented here, these models utilize techniques to reduce, 
and therefore, simplify a higher order system. While 
each of the papers referenced present good simulation 
and implemented results compared to actual data, in our 
work we are concerned with validating that the reduced 
order observer performs better or just as well as the full 
order observer. Moreover, for ICEs specifically, 
reducing the size of the observer without reducing the 
size of the model beforehand has never been investigated. 
There has been much work in recent years to model the 
dynamics of the internal combustion engine (ICE) [18- 
23]. The ICE has become so complex that the engine is 
usually split into subsystems, each with its different tasks 
and constraints. A critical aspect of these tasks related to 
state estimation can provide the best fuel economy and 
emissions efficiency. Another important task is finding 
ways to reduce manufacturing cost by reducing the 
amount of redundant components. Such components 
include sensors, which can be unreliable, expensive and 
sensitive to environmental changes (pressure, 
temperature, water etc...). Indeed, electronic control 
units (ECUs) used in modern powertrain systems in 
today’s vehicles already depend on well-designed 
observers to estimate states to work with state feedback 
control. 

While the reduced order observer’s ability to estimate 
and to control is well documented, it has not been utilized 
in engine models. In this paper, a full order observer is 
designed to estimate all the states of a third order ICE 
model. A comparison is then made between the 
performances of the full order observer to that of reduced 
observers. A simple step input command is used for all 
observer simulations and we analyse the steady stead 
behaviour of the responses. Through comparison of 
these observers’ performance, we validate that a reduced 
order observer can perform better or just as well as a full 
order observer. 

In section 2, an outline of the reduced observer is given. 
Section 3 provides the design and implementation of the 
reduced observer on the ICE model. In section 4, an 
analysis of the results is discussed and corresponding 
simulations are provided. Section 5 is a conclusion of the 
work done in this paper. 
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2. Reduced Order Observer Theory 

The design of a full order observer is well known and can 
be found in [1]. Here we outline the design for the 
reduced order observer used. 

For a Linear Time Invariant (LTI) system written as: 



Ax(t ) + Bu(t ) 


(1) 


Cx{t) 


(2) 



derived. 






w = 


Tx 


(7) 


x = 


T' l x 


(8) 


Tx = 


TAx + TBu 


(9) 


w = 


TAT~ l w + TBu 


(10) 



where A, B and C are the state, input and output constant 
matrices, respectively. A, B and C have dimension 
fl xfl, r xTl and TTl x U, respectively, H is the number 

of states in the system, r is the number of inputs and 
Hi is the number of directly measurable states used as 
outputs. A transformation matrix, T, is chosen, 



T 



C 

D 



(3) 



where D is a matrix that needs to be chosen to make the 
transformation matrix the appropriate dimensions and 
non- singular. Using prior knowledge of the structure of 
the output matrix, C, matrix D is chosen and the 
transformation matrix can be written as: 



The LTI system in (1-2) is multiplied by the 
transformation matrix, T, to get (9). Manipulation of (7-9) 
creates the new transformed LTI state equation: 

W = Aw + Bu (11) 

where A and B are the new state and input matrices, 

respectively. The new state equation can be written as 
follows, 





j 
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w = 


= A 




V 




V 



( 12 ) 



The vector v contains the unknown states, which need to 
be estimated. It is advantageous at this point to partition, 

^ /V 

A and B of the transformed state vector as shown here: 



~C~ 
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^ n-m 
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_o m 


^ n-m _ 



(4) 



with I m and I n _ m being identity matrices indexed by 

their respective dimensions. Similarly, 0 n _ m and 0 m are 
zero matrices indexed by their respective dimensions. 
Matrix D is chosen to be the ( Jl — m)xin matrix 

KL / , which ensures that the transformation 

matrix is always non-singular, as long as the output 
matrix, C, is the TTlxTl matrix [/ m 0 /? m ] . A new 

state vector is then derived using the transformation 
matrix, 



w = Tx 



(5) 
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A 2 
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(13) 



Ajp Aj 2 , A 21 and A 22 are matrices with dimensions 
mxM , mx(n-m ), ( n-m)xm and 

(n-m) x(n-m) , respectively. B x and B 2 have 

dimensions Hi x r and ( H — Hi) x r, respectively. The 
expansion of (13) gives two separate equations: 



y = A n y + A 12 v + B x u (14) 

/v zv zv 

v = A ll y + A 22 v + B 2 u (15) 



Substituting (3), the new state vector can be written as a 
relationship between the system output, y , containing 

m measurable states and a new matrix with ( Tl — Hi) 
states, 



w = T x = 


"c" 


X = 


J 




D 




V 



Vector v contains the remaining state variables, which 
need to be estimated with the use of an observer. This 
observer has reduced dimensions because the measurable 
states from the state vector are not included. The 
observer will have the form: 

v’ = (y 2 -GA l2 )v +G(y-A n y-B 1 u) 

(16) 

+•^21 y + B 2 u 



V is a (n — fri) x Tl matrix containing the states which 
need to be estimated. Now that the new state vector is 
derived, a new transformed linear time invariant state 
equation is derived. By manipulating (5), (7-8) are 



The following relationship, Z — V — Gy , is used to get 
rid of the y term in (16). G is the observer gain of 
dimensions ( H — Hi) x Hi. The gain can be found using 
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gain determination methods such, as a Linear Quadratic 
Regulator (LQR), which calculates a gain that will 
minimize an error criterion as discussed in [24-25]. 
Another method to determine the gain is pole placement, 
which places the observer poles in the stable region of the 
s-plane to ensure the observer estimation converges to the 
state [26-27]. For this paper, we used the pole placement 
method. Even though the LQR method is more optimal 
since it minimizes an error criterion, the pole placement 
method is simpler in practice. It can be used to achieve 
control by simple trial and error and using certain rule of 
thumbs. However, the pole placement method has its 
limitations and its lack of ability to optimize the gain can 

lead to low observer performance. V and V indicate the 
estimated state variable vector of v. Getting rid of the 
differentiator term y gives the final observer form: 



Z 



( A 22 -GA n )z + ( A 22 -GA l2 )Gy 
+ (A 21 - GA n ) y + (B 2 - GB X ) u 







"-35 
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0 " 


* 




" 35 


*2 
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-3.822 


-6.552 


*2 
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A- 




0 


3.2 


1.84_ 


— * 3 - 




_ 0 



y 



[i i i] 




( 19 ) 



In order to present the fundamental stability of the model, 
we calculate the controllability and observability of the 
system. It is known in literature that by taking the 
absolute value of the determinant of the respective 
matrices, we are able to get a measure of just how 
controllable and observable the system is. The results are 
presented in Table 1. The absolute values of both 
determinants, presented in columns 2 and 4 of Table 1, 
yield two large numbers that are far from 0. This 
represents a system that is highly controllable and 
observable. 



The reduced order observer model in (17) is now 
obtained as a linear relationship between the measurable 
states, y , the measurable input, U , and the estimated 
state vector, z , which contains the remainder of the 
system’s states. 




Figure 1 . Reduced Order Observer Block Diagram 

Figure 1 depicts a block diagram of the reduced order 
observer. 

3. Design and implementation 

The third order ICE model along with its parameters 
were created and thoroughly outlined by Baumgarther 
and Geering in [28]. The model, of a Chrysler Neon, was 
slightly adjusted for this paper. The inputs, u , were 
reduced from three to one. The state vector is as follows: 

x x : Throttle Angle ( Degrees ) 

x 2 : Intake Manifold Pressure ( Torr ) 

jc 3 : Engine Speed ( RPM ) 

u : Commanded Throttle Angle ( Degrees ) 

The mathematical model for nominal parameters is 
presented as: 



Observability Matrix (N) 


ldet[N]l 


N = 


o o r 

0 3.2 1.84 

13.478 -6.34 -17.581 




-43.131 


Observability Matrix (M) 


ldet[M]l 


M = 


"35 -122.5 42.87" 

0 147 -57.23 

0 0 472 _ 
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Table 1. Observability and controllability of ICE model 



Full Order Observer design: First, a full order observer 
is designed to estimate all the states of the system. A 
standard Luenberger observer is designed using the form: 

x(t) = Ax(t) + Bu(t) + G x [y(0 - K0] (20) 

m = Cx{t) (21) 

The parameters for (20-21) are presented in Table 2. x 
denotes the estimated state vector. The full observer gain 
which we shall refer to as Gi was found using the pole 
placement method and placing the full observer poles (A 
- GiC) at a location that would make them faster than the 
system’s poles. To ensure the observer stabilizes well 
before the system, we made the observer poles three 
times faster. The corresponding observer gain and 
system’s poles are also presented in Table 2. 

Reduced Order Observer design (2 nd Order): Next, we 
look at the case where the throttle command can directly 
be measured so there is no need to estimate it using an 
observer. The reduced observer, then, becomes a second 
order as opposed third order as the full observer. The 
transformation matrix (4) was chosen along with the 
partitioned state and input matrices. Next, we find the 
observer gain, G 2 , using the pole placement method. 
Similarly to the full order observer, the poles of the 
second order reduced observer (17), which 
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are(A 2 2 — G 2 A 12 ) , were places at a location that makes 

them three times faster than the system’s poles. All of 
the aforementioned parameters are presented in Table 3. 



State Observer 
model system 

parameters 


Parameter matrix 


State Matrix 


A = 


”-35 0 0 ' 

4.212 -3.822 -6.552 

0 3.2 1.84 




Input Matrix 


B = 


35 

0 

0 




Output Matrix 


C = [ ill] 


Gain Matrix 


G l ~ 


"2.1133 " 
3.5262 
26.7216 




System Poles 


A = 


-35 1 

_-0.991 ± 3.598 lj 



Table 2. Full observer design parameters 
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model system 
parameters 
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Gain Matrix 
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Transformation 

Matrix 


T = 


T o o' 
0 1 0 
0 0 1 





Table 3. 2 nd Order observer design parameters 



Reduced Order Observer design (First Order): Lastly, 
we look at the case where both the throttle command and 
the manifold air pressure are both directly measurable. 
There is now only a need to estimate the engine speed. 
The reduced observer is now in its simplest form as a 
first order observer. The transformation matrix is the 
same as for the 2 nd order reduced observer. The new 
partitioned state and input matrices have different 
dimensions as the 2 nd order reduced observer. The 
observer gain, G 3 , is again found using the pole 

placement method. The poles of (A^ — G 3 A 12 ) were 

placed at a location that makes them three times faster 
than the system’s poles. All of the parameters are 
presented in Table 4. 
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Table 4. 1st Order observer design parameters 
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4. Results and Analysis 

In this section, the simulation results are presented and 
analysed. First, we take a look at the simulations of all 3 
observers; Full order, 2 nd order reduced and 1 st order 
reduced. Figure 2 depicts all three states, which were 
estimated using the full order observer along with the 
actual measured values of the system. Figure 3 shows 
the result of the second order reduced observer states; 
manifold pressure and engine speed, and compares it to 
their corresponding full observer estimate and actual 
measurements from the systems. Figure 4 presents the 
simulation results of the first order reduced observer 
where the only state is the engine speed. The full 
observer and second observer estimate of the engine 
speed are also included for comparison, along with the 
actual measured speed from the system. 



r 



"Actual Throttle Angle 
Full Observer Throttle Angle 



Time (s) 




Figures 5-7 depict the error dynamics of each of the 
observers: full, 2 nd order reduced observer and 1 st order 
reduced observer. 



(0 20 



Full Observer Error Dynamics 
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Figure 5. Error dynamics of full observer: Commanded throttle position, 
Intake manifold pressure and Engine speed 



Reduced Order Observer (2nd Order) Error Dynamics 
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Figure 3. Simulation results of full observer and 2 nd order reduced 
observer 




Figure 4. Simulation results of full observer, 2nd order reduced 
observer and 1 st order reduced observer 
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Figure 6. Error dynamics of 2 nd order reduced observer: Intake manifold 
pressure and Engine speed 



Reduced Order Observer (1st Order) Emir Dynamics 




Figure 7. Error dynamics of 1 st order reduced observer: Engine Speed 

We take a special interest at the error dynamics to 
evaluate the performance of the reduced observers. All 
observers have been very well designed and visually their 
error converves towards zero asymptotically and quickly. 
Comparing the error dynamics for the manifold pressure 
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in the full and 2 nd order reduced observer, it is evident 
that the reduced observer converges faster than its 
counterpart. It does seem , however, that the 2 nd order 
reduced observer’s average error thorughout the entire 
simulation is higher than that of the full observer. In the 
same fashion, we take a look at the error dynamics of all 
three observers for the engine speed estimation. Again, 
in all three observers, the estimation converges quickly 
with the 2 nd order reduced observer having the least 
initial error, and the full observer converged fastest. The 
1 st order reduced observer, nonetheless, was only about 
250 ms slower to converge which is an insignificant 
amount of time. 



5. Conclusion 

In this paper, we addressed an issue related to state 
estimation of typical internal combustion engines. This 
issue involves the redundancy, high costs and 
unreliability that sensors are subject to and can cause. 
We validated the use of reduced observers as a solution. 
A reduced observer can estimate the unmeasurable states 
using measurements of those which are measured directly. 
This reduces the order of the observer, which otherwise 
would need to estimate more states. In practice, this 
allows for the removal of sensors. A third order ICE 
model was used for which we designed and implemented 
a full order observer and two reduced order observers to 
estimate its three states. The two reduced order observers 
were of second and first order, for the cases where 1 and 
then 2 states are directly measurable. An analysis of the 
error dynamics of all three observers was discussed. The 
reduced observers converged just as fast as or even faster 
than the full observer. Also, the reduced observer 
showed an insignificant amount of steady state error 
which was only slightly larger than or just as low as the 
full order observer. This demonstrates that the use of 
reduced order observers can be extended further to 
present day powertrain control modules, which already 
use observers in their algorithm. They can be utilized to 
cheapen cost, reduce unreliability and decrease system 
complexity. 
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Abstract 

The paper discusses the Simulink model of lithium 
polymer battery. Model parameters are temperature and 
state of charge dependent. This model is two branch 
(parallel R-C) model. 2-D Look-up table method is used 
to represent the relation between model branch 
parameters and state of charge. Matlab user defined 
function is used to calculate the open circuit voltage and 
internal dc resistance. One whole block is dedicated to 
state of charge calculation. This block is able to 
autonomously detect the transition from charge to 
discharge and vice versa. Also if charging or discharging 
current stops flowing, it resets the current integrator 
initial condition to the last state of charge calculated. 
Comparison between experimental and model results 
showed about 3% maximum error in dynamic (interval 
test) response and 5% error in static response of battery 
discharging (until 10% state of charge is left in the 
battery). 

Keywords: Battery, look-up table, state of charge. 

Nomenclature 

PSCAD Power System Computer Aided Design 

R-C Resistor- Capacitor 

SOC State Of Charge 

DAQ Data Aquisition 

A/D Analog to digital 

D/A Digital to analog 

AH Amp-Hour 

V t , v Terminal voltage 

V oc , E Open circuit voltage 

R dc , R 0 Internal DC resistance 

Ri, R s Resistance of short time constant 

R 2 , Ri Resistance of long time constant 

Ci, C s Capacitance of short time constant 

C 2 , Ci Capacitance of long time constant 

I Current 

t Time 

T Constant time interval 

k Constant of temperature directly proportional to 

temperature (Table-5) 

a m Coefficients of quadratic equation (Table-6) 

Q Battery capacity in amp- sec (AS) 



1 Introduction 

Increasing the throughput of electric vehicles (EV) is 
very important. To do that, main component of the EV, 
the battery, must be tested, analyzed and modeled before 
it can be installed in EV. Dynamic models of batteries are 
important for analysis, design and performance of EV [1- 
3]. Most models developed until now are for small 
batteries (low capacity), used in cell phones, computers, 
and portable electrical appliances which show less 
nonlinearity. But electric vehicle grade batteries (high 
capacity) are highly nonlinear in nature. The surface 
kinetics, thermal behavior, reaction rate of big batteries 
differs from small batteries. 

Proposed model is adopted from [4] (which is developed 
based on the experiments on batteries used for laptop 
computers, PDAs, mobiles etc.) and it is based on 
equations, so it can be implemented in different 
simulation platforms as it has been implemented in 
PSCAD [5-6]. From all the models [1-23], model from [4] 
showed the best accuracy in dynamic as well as static 
characteristics. For the development of model presented 
here, each branch voltage (voltage across dc resistance, 
voltage across each pair of R-C, open circuit voltage) 
calculation is divided into blocks. Developing a model in 
different blocks is to identify different function of each 
block that develops the model. Each block represents 
contribution to terminal voltage calculations. These block 
parameters depends on state of charge (SOC) and 
temperature. The parameters are polynomial equations 
depending on SOC and temperature [7]. 

Wide variety of models has been developed by 
researchers with varying degree of complexity. Electro- 
chemical models developed in [8-10] mainly designed for 
optimization of physical design of the battery. But these 
models involve complex time-variant partial differential 
equations which require days to simulate and requires 
complex algorithms. Mathematical models [11-13] lacks 
in practical meaning but useful to system designers, uses 
mathematical methods like stochastic approaches [12] to 
predict battery runtime, efficiency or capacity. However, 
mathematical models do not give out the SOC-V 
information that is important in circuit simulation and 
optimization. Moreover, mathematical models provide 
inaccurate results in the order of 5%-20% error. 
Electrical models [1-7, 14-23] provides the accuracy 




33 





ACSE Journal, Volume 14, Issue 1, ISSN 1687-4811, ICGST LLC, Delaware, USA, August 2014 



somewhere in between electro-chemical models and 
mathematical models (around l%-5% error). Following 
few paragraphs describes few of the developed electrical 
models. 

In figure (1), a Thevenin’s equivalent circuit model is 
presented, which has been extensively used to model the 
discharging dynamics of the batteries such as Lead acid, 
Lithium polymer, Lithium -Ion (Li-Ion) etc. [1 -3,5-7]. 
The easiest model is the first order Thenvenin’s model 
with one parallel R-C branch and its parameters can be 
easily estimated experimentally. The R-C combinations 
predict the response to a transient load at specific SOC’s 
by assuming that the open circuit voltage is constant. 
Therefore, this model is unable to reflect the influence of 
the SOC to the battery behavior properly. Taking this into 
consideration the battery voltage response might be 
different from first order model response [14]. So, higher 
order model is required to capture the behavior of the 
battery. For some batteries second order model is 
sufficient e.g lithium polymer batteries. 



— VvV 





Figure 1 : Thevenin’ s equivalent model for discharging 
battery [14] 



An impedance based circuit acquires an AC-equivalent 
impedance model in the frequency domain. It then uses a 
complex equivalent system (Z ac ) to fit the impedance 
ranges. This fitting process is very difficult and complex. 
Additionally, since these models only work for fixed 
SOCs and temperatures, they cannot predict battery 
runtime or DC response [15]. 

Runtime based models use complex circuits to simulate 
the DC voltage response and runtime of the battery. This 
model used two different circuits to depict the battery 
characteristics. One circuit depicts a capacitor, which 
contains the value of the battery capacitance, and a 
dependent current source which represents the batteries 
state of charge. And the second circuit consists of the R- 
C networks which represent the relationship between the 
battery current and terminal voltage [4]. This circuit is 
similar to the Thevenin’s equivalent circuit. But these 
models do not take into consideration the changing 
behavior of passive elements with varying SOC. 
Developed model incorporates the effect of SOC and 
temperature on different parameters of the battery model. 
Also, the model is valid for different chemistries 
provided the equations for different parameters of the 
battery depending on SOC and temperature are available, 
but at the time model is tested for Li-Ion chemistry. 

The remainder of the paper is organized as follows: 
section (2) describes the battery evaluation lab, section (3) 
focuses on test system and test method to extract the 
parameters of the battery model, section (4) emphasizes 
on the development of the model in 

MATLAB/SIMULINK, section (5) emphasizes on SOC 
calculation block, section (6) analyses and compares the 



results obtained from the simulation to the experimental 
results and section (7) concludes the paper with 
emphasizing the contribution of the paper in developing 
the battery model with parameters depending on SOC 
and temperature. 

2 Battery Evaluation Lab 

The tests on the batteries are carried out at the battery 
evaluation lab at University of Massachusetts Lowell. 
The lab is equipped with current regulators controlled by 
computer. There are two types of regulators. 

1. 8 current regulators with capability of upto 10V 
and 320 amps source or sink. 

2. 4 current regulators with capability of upto 20V 
and 160 amps source or sink. 

All the current regulators are capable of going from 
charging to discharging instantaneously [16] and are 
controlled via A/D and D/A controllers, which gets the 
commands from data acquisition (DAQ) system 
controlled by computer. All the programs are written in 
BASIC. Different current profiles can be easily written 
with BASIC programming language. You can load 
external current command file, so that realistic load tests 
can be done. 

Lab is also equipped with two environmental chambers to 
test the batteries at different temperatures. Range of 
temperature control is from -20°C to +50°C. Before every 
test the battery is soaked to ambient temperature. If the 
temperature is 0°C or below battery is soaked overnight 
to make sure that center of the battery is at desired 
temperature. 

3 Test system and test method 

Test system is shown in figure (2). As it is seen, the 
battery which resides in the environmental chamber is 
connected to the current regulators and MUX. Current 
regulators and environmental chamber are controlled 




Figure 2: Test setup 



The battery model must resemble the dynamic and static 
characteristic of battery behavior in real world. So, two 
types of tests are done to extract the parameters of battery 
model. Interval test is done to monitor the battery 
dynamics and provides information about the short time 
and long time constant parameters for the battery model. 
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Constant discharge test is done to obtain the DC 
resistance. 

Interval test is constant discharge followed by rest. 
Constant discharge in interval test is pre specified and set 
to lOAh (0.1 SOC). The rest period is set to 10 minutes 
so that by that period the dynamic effect of short time 
and long time constants is stabilized and battery reaches 
to open circuit voltage. Figure (3) shows the interval 
discharge at 0°C and 100 amps. For constant discharge 
test, battery is subjected to constant current of lOOamps 
at different temperature. Figure (4) shows the constant 
discharge at 20°C and 100 amps. 

Presented model is developed for discharging battery, but 
the same tests (interval test and constant discharge test, 
explained in detail later in the paper) can be done to get 
the parameters of the battery during charging and can be 
incorporated in the same model. Same parameters are 
used for charge and discharge for this model. 




Figure 3: interval discharge test done at 0°C and 100 
amps 




Figure 4: constant discharge test done at 20°C and 100 
amps 



To be able to predict the behavior in different weather 
condition tests are done at different temperature settings. 
The temperature range is from 0° to 40° C. 

Parameter extraction: Thevenin second order circuit 
model solution is given by: 

t_ 

v t = V oc - % - IRi (1- e RA ) - IR, (1- e^) (l ' 

The method used to extract the parameters for model is 
adopted from [14]. It’s a simple analytical method. In this 
paper explicit formulas are used to compute the nth-order 
model. But we used its 2 nd order model derivation. We 
tweaked the results of [14] to work for interval test. 
Instead of being able to work at particular SOC 
beginning, we added a loop that at any instance of time if 



current goes to zero and after some time discharging 
current starts to flow again program recalculates the 
parameters with new values of (v, t) specified at constant 
time intervals T. Tables 1-4 show the parameters 
extracted at different temperature. 

Figure (5) and figure (6) shows the graphs of open circuit 
voltage and DC resistance at different temperature 
respectively with both the measurement and simulation 
results. 



Table 1 R s -SOC values at different temperatures 



SOC 


Temperature (°C) 


0 


20 


40 


0 


0.00086 


0.0006 


0.00016 


0.1 


0.000722 


0.00025 


0.000122 


0.2 


0.000621 


0.00022 


0.000121 


0.3 


0.00062 


0.00022 


0.00012 


0.4 


0.00050061 


0.00020061 


0.00010061 


0.5 


0.000222291 


0.000162291 


0.000102291 


0.6 


0.00020408 


0.00013408 


0.00010408 


0.7 


0.00018078 


0.00010078 


0.00010078 


0.8 


0.00017836 


0.00025836 


0.00011836 


0.9 


0.000203084 


0.000063084 


0.000103084 


1 


0.00014976 


0.00012976 


0.00008976 



Table 2 C s -SOC values at different temperatures 



SOC 


Temperature (°C) 


0 


20 


40 


0 


20000 


56000 


60000 


0.1 


22011 


40011 


60011 


0.2 


26379 


40379 


60379 


0.3 


31000 


71000 


61000 


0.4 


51696 


81696 


61696 


0.5 


57680 


107680 


67680 


0.6 


54089 


44089 


64089 


0.7 


55760 


225760 


65760 


0.8 


57760 


187760 


67760 


0.9 


56265 


66265 


66265 


1 


52847 


202847 


62847 



Table 3 R r SOC values at different temperatures 



SOC 


Temperature (°C) 


0 


20 


40 


0 


0.00061 


0.02 


0.00026 


0.1 


0.00057 


0.02 


0.00025 


0.2 


0.00053 


0.00044 


0.00023 


0.3 


0.00049 


0.00047 


0.00022 


0.4 


0.00047 


0.00047 


0.00027 


0.5 


0.00042 


0.00045 


0.00025 


0.6 


0.0003 


0.00039 


0.0001 


0.7 


0.00016 


0.00023 


0.0001 


0.8 


0.00017 


0.00007 


0.00017 


0.9 


0.00015 


0.000155 


0.00017 


1 


0.00016 


0.00017 


0.0001 
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Table 4 Q-SOC values at different temperatures 



SOC 


Temperature (°C) 


0 


20 


40 


0 


249350 


249350 


249350 


0.1 


314350 


164350 


314350 


0.2 


329880 


419880 


419880 


0.3 


352000 


402000 


402000 


0.4 


401950 


401950 


401950 



Open circuit voltage for different temperature 

4 . 25 1 1 1 1 1 1 1 1 1 




State of charge (SOC) 



Figure 5: Open circuit voltage of the battery at different 
temperature with measurement and simulation 



0.5 


401240 


401240 


401240 


0.6 


401660 


401660 


401660 


0.7 


401870 


401870 


401870 


0.8 


408170 


408170 


408170 


0.9 


409930 


409930 


409930 


1 


408920 


508920 


408920 



In the interval discharge test open circuit voltage (V oc ) 
and internal dc resistance (R dc ) is recorded at each end of 
rest period: 



K c = v(0) (2) R dc = v(0)-v(0 + )/I (3) 

After finding V oc and R dc at different SOC points, 
polynomial equation is derived using curve fitting 
toolbox to relate V oc and R dc with SOC. With three 
different temperatures we have three different equations 
for V oc and three different equations for R dc . Combining 
these equations produced following equations: 



7 (4) 8 (5) 

V oc =E b n S0C " \ = X>„soc" 

n=0 n=0 




Figure 6: DC resistance of the battery at different 
temperature with measurement and simulation 

Table 5 Relation between temperature and constant k 



Temperature (°C) 


k 


0 


0 


20 


2 


40 


4 



Where, 

b„ = £> m k m )} n 

m = 0 

k and a m values are presented in table 5 and table 6 
respectively 

At 0°C the dc resistance rises very quickly toward the end 
of discharge. Also values of dc resistances at different 
SOCs are relatively higher than 20°C and 40°C, which 
happens to be due to slow kinetics at low temperatures. 
As temperature increases the value of dc resistance 
decreases. 

4 Simulation model 

The model is implemented in MATLAB/SIMULINK 
simulation environment. Visualization of equations 
becomes so easy with Simulink’s graphical presentation. 
User defined function written in Matlab script editor can 
be used in SIMULINK. 

This model is two circuit model connected through 
nonlinear voltage controlled voltage source and a current 
controlled current source. The first circuit consists of big 
capacitor representing SOC of the battery and the other 
circuit is two pairs of parallel R-C branch connected in 
series with internal dc resistance [17] as shown in figure 
(V). 



Battery Lifetime Voltage-Current Characteristics 

i 1 t 1 




Figure 7 : Proposed battery model adopted from [4] 



Self-discharge resistance is not considered in the circuit 
as the battery is subjected to frequent use. R dc and V oc can 
be found from constant discharge test, but it is only 
possible for initial condition e.g if the constant discharge 
test starts at 90% SOC than R dc and V oc at 90% SOC is 
available. But to establish the relationship of R dc and V oc 
with SOC, interval test is done. From this test other 
component values are obtained such as R s - C s and R r Q 
components of short time constant and long time constant 
respectively. Battery behaves differently under different 
climates, so relationship of R dc , V oc , R s , C s , Ri, Q with 
SOC changes at different temperatures. All the model 
components and SOC are related by polynomial equation. 
So while the model is running, model calculates the 
values of all the components dynamically and uses it to 
find the terminal voltage at that particular SOC. So 
overall all the calculations depend on SOC of the battery. 
One block is allocated for the calculation of SOC, which 
takes care of following actions: 

1 . Sudden start of charge from discharge and vice- 
versa 

2. Rest after charge or discharge. 

During these conditions the integrator resets itself to 
new initial condition. 
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Table 6 Values of coefficients a 2 , ai, a 0 of different b, which are coefficients of Y oc and R dc 





Voc 


Rdc 


a 2 


ai 


a 0 


a 2 


ai 


a 0 


b 8 


- 


- 


- 


0.01194875 


-0.1033125 


0 


b 7 


7.3705 


-14.741 


0 


-0.05989975 


0.4332045 


-0.075432 


b 6 


-29.075 


77.1571 


-42.848 


0.10599375 


-0.6961025 


0.22984 


b 5 


48.6737 


-167.355 


153.79 


-0.07906375 


0.5003425 


-0.20603 


b 4 


-45.1448 


192.389 


-219.71 


0.017924363 


-0.1190662 


0.0082249 


b 3 


25.0809 


-125.036 


159.4 


0.00658875 


-0.035942 


0.075254 


b 2 


-8.2449 


45.3 


-61.088 


-0.00388465 


0.02335445 


-0.035913 


bi 


1.4449 


-8.3407 


11.983 


0.00077195 


-0.0046413 


0.0083023 



All the equations for the model of discharging batteries 
are as follows: 

T n 1 

SOC(t) = SOC(O) - J — xl(t)dt 
o Q 

Where, 

SO C(0) -initial SOC 



V, = V oc - - n\ (l-e R,c ’)- n^(l-e"^) (7 ' 

In the figure (8), source (supply) named “current” can be 
set to run in constant discharge mode or interval 
discharge mode. Charging can also be incorporated in the 
same source. This current is fed to four subsystems 
namely SOC (scion colored block), Internal DC 
resistance (grey colored block), short time constant 
(green block), and long time constant (orange block). 
First the SOC block is executed. Value of temperature is 
fed to 1-D look up table called k and also other 2-D look 
up tables R s , C s , R\ and Q. Using SOC and temperature 



values of R s , C s , Ri and Q are determined, and fed to 
short and long time constant blocks. V oc and R dc are 
written in user defined matlab function. Values of k and 
SOC are supplied to this function to evaluate the 
equations. Vt is determined in four parts. The four parts 
are as follows: 

1 . Open circuit voltage V oc 

2. Internal DC resistance block is the mathematical 
representation of IR dc 

3. Short time constant block is the mathematical 



4. 



t_ 

representation of ^ (l—e *&) 

long time constant block is the mathematical 
representation of ^^“4) 



To get the value of V t parts 2, 3 and 4 are subtracted 
from 1. Two conditions are applied if the voltage goes 
above 4.2V or voltage goes below 2.7V the model stops 
execution. 




Figure 8: Simulation model of lithium battery 
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5 SOC calculation 




There are three modes of operation during which SOC 
varies accordingly, 

1. charge mode, SOC increases during this mode 

2. discharge mode, SOC decreases during this 
mode 

3. rest mode, SOC remains unchanged during this 
mode 

SOC calculation is a big part of this model and it is 
done in SOC block (figure (9)) of the model. Battery 
capacity, initial SOC and current are input to this model 
and SOC is the output. First the capacity is converted into 
amp- sec (AS) format from amp-hour (AH) and inverted 
to be multiplied by current. Reset value of SOC is added 
to compensate for full discharge or charge condition as 
well as transition from one mode to another. So if battery 
is discharging and SOC reaches to 0 the output of the 
adder is set to 0. And if the battery is charging and SOC 
reaches to 1 the output of the adder is set to 1. Reset 
value of SOC comes from block named “corrected SOC 
after full discharge or charge”. If the battery is in 
discharge mode and it reaches to 0 SOC but the battery 
does not go in charge or rest mode than this block holds 
the 0 SOC, same goes for charge process. But this block 
does not have to hold these conditions for long because at 
0 SOC output voltage is 2.7V and at 1 SOC output 
voltage is 4.2V so, simulation stops. Also this block 
works as memory i.e, if battery goes into rest mode or 
charging mode from discharge mode, this block 
memorizes the value of SOC at the time of transition and 
uses it as initial SOC for next mode. After the addition 
block one more condition is checked if the current is zero 
than the same value of SOC is maintained. 




6 Results 

For the interval test current value has to be equal for all 
temperatures so that we don’t have to change the time 
value for discharge during each interval. Figure (10) 
shows the current profile for interval discharge test. 
Positive current is for discharge and negative current 
represents charging process. Interval test was done at 
different temperatures. 




Time (sec) 



Figure 10: current regime for interval discharge test 

The current profile is chosen such that battery discharges 
for 10% SOC and then enters into rest mode. For tested 
battery of 100 Ah 10% SOC corresponds to lOAh. So if 
the battery is discharged at 100 amps current then the 
interval should be of 360secs to discharge the battery 
lOAh. Figure (11 (a)) and figure (11(b)) shows the 
terminal voltage comparison between measurement and 
simulation of interval tests at 0°C and 40°C, lOOamps 
respectively. Figure (12) shows the terminal voltage 
results of measurement and simulation of constant 
current discharge tests at 0°C and 40°C, lOOamps 
respectively. 
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Terminal voltage @ 0 deg celcious 
and 100 amps discharge rate 



Terminal voltage @ 40 deg celcious 
and 100 amp discharge rate 



<D 

CD 

2 

o 

> 




Time (sec) 
(a) 



Time (sec) 

(b) 




Figure 11: Terminal voltage of the battery subjected to transients of the interval of 10% SOC followed by 10 min rest (a) 
at 0°C and 100 amps discharge rate (b) at 40°C and 100 amps discharge rate. 

components of short term time constant and long term 
time constant. The polynomial equations of open circuit 
voltage and DC resistance are of the 7 th and 8 th order 
respectively to capture the precise value. This model is 
designed for both charge and discharge. This model 
makes use of tabular method to find values of 
components of short term and long term time constants 
depending on SOC. Interpolation and extrapolation 
capabilities of the table used in Simulink help in finding 
values other than test values. 

Results show that the maximum error between model and 
experimental data between 100%- 20% SOC is about 
<1% for interval test and 1.2% for constant discharge test. 
But at the end of simulation (20%-0% SOC) this error 
increases. Because the interval test and constant 
discharge tests are done at different open circuit voltage, 
error in the terminal voltage is high for constant 
discharge test. At lower temperatures dc resistance 
increases unexpectedly at low SOC. 
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Abstract 

This paper presents a new hybrid dynamic approach to 
describe and control sequential biped robot motions 
including sitting and standing. Transitions between these 
motions are expressed as a hybrid dynamic system that 
express each pose of the biped. Control of the biped in 
each state is performed using nonlinear backstepping 
with integral action. Stability and robustness of the 
ensuing biped robot control are addressed. Simulations 
are used to further demonstrate the effectiveness of 
control. 

Keywords: nonlinear control, integrator backstepping, 
hybrid dynamic systems, biped robot. 

Nomenclature 

q t Discrete states of a hybrid system 
X Continuous-time state space of a hybrid system 

D Domain of continuous states in a discrete state q 

f Function describing the continuous dynamics in a 
state q . 

E Edges, possible transitions from states 
G tj Guards, conditions that must be met to transition 
to another state 

R Resets, reset maps between states 
D (q ) Inertia matrix of a robot, size of nxn. 

H [q , q ) Friction, coriolis and centripetal forces, size of 
nxl. 

G [q ) Gravitational forces acting on each link’ s center 
of gravity, size of nxl. 

u Actuator torques acting on each link, size of nxl. 
S y HDS describing a sitting biped 

^2 HDS describing a standing biped 

c l’ "' ,c 2 n Backstepping controller parameters for an nxn 
robot. 

k ^, ..., k n Backstepping controller integral gains for an 
nxn robot. 



1. Introduction 

There are several aspects to walking bipedal motion 
including the multi degree of freedom dynamic model of 
the biped [1], control and coordination of the walking 
gait [2]-[4], actuator control, and disturbance rejection 
because of terrain and other external forces. Each of 
these topics is very important by itself and combining 
them together to form a working system is a continuing 
challenge. 

In this paper, a framework for hybrid dynamic systems 
(HDS) is presented to show that it is desirable to organize 
biped robotic motions into multiple levels of subsystems 
for control. Organizing biped motions into hybrid 
dynamic systems, sometimes with a reduced model order, 
allows for greater understanding when controlling 
transitions such as sitting and standing or between 
standing and sitting. Examples are shown by using 1-link 
and 2-link representations of a sitting biped robot model 
derived using Lagrangian equations of motion [5]. 

Around the hybrid states the biped is controlled using 
nonlinear backstepping control. This method of control 
has several benefits including robustness against 
unmodeled dynamics and meeting stability and 
performance requirements [6]- [12]. A scheme for 2-link 
robotic actuators will be applied to the biped system. 
Each HDS has a unique controller and the transitions 
between HDS systems require the controller’s smooth 
transfer under predetermined conditions. Sequencing of 
HDS states and control is shown through examples of 
motions from a sitting biped to a standing biped. 
Implementations of the simulations of models are 
performed in Matlab/Simulink. Simmechanics is used to 
model the biped robot. Zhao and Liu proposed a 3D 
Simmechanics model of a biped for simulation which 
was capable of using joint angle actuator drivers to 
simulate the kinematics of the biped [13]. However, their 
model is not able to simulate dynamic feedback control. 
In this paper a simplified 2D Simmechanics model is 
presented that is capable of simulation with actuator 
torque for dynamic feedback control. 

While hybrid dynamic systems and backstepping use is 
well documented in many applications, both have not 
been combined and used extensively in biped robotics. 
In this paper a HDS system and transitions are designed 
to describe biped robot motions. In section 2, an example 
of HDS definitions for sitting and standing motions is 
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given. Section 3 provides the biped model derivation 
using Lagrange Equations. In section 4, backstepping 
controllers with integral action are designed for sitting 
and standing states. Section 5 shows simulations of the 
sitting to standing events, compares the backstepping 
control to a classical PD control [14] and explores 
stability and robustness. Section 6 is a conclusion of the 
work of this paper. 

2. HDS Model Of Action Motions 

Hybrid Dynamic Systems contain both continuous 
dynamics and overlaid discrete events. Within a hybrid 
state, the continuous dynamics are described by a system 
of differential or difference equations. While in a hybrid 
state the dynamics flow in time but can transition to the 
same state or to another state when certain guard 
conditions or constraints are met. Jumps to other hybrid 
states may cause state variables to reset based on defined 
jump event maps. The order of continuous dynamics in 
each hybrid state may also vary and the maps will define 
how this change and its effects are handled. This 
notation of an HDS system is known as the hybrid 

automaton of the form S = {Q,X ,D, / , as 

shown in Figure 1 [15]- [17]. 



G n (q l ,q 2 ) x + eR n (q l ,q 2 ) 




Figure 1. Hybrid Dynamic System model with two hybrid states. 

The variables in Figure 1 are defined as follows. 

q t - Discrete states of the hybrid system (mode of the 

system); e.g. falling, rising, accelerating, decelerating. 

X - Continuous-time state space of the hybrid system. 

D - Allowable domain of continuous states in a 
discrete state q . 

/ - Vector field function describing the continuous 

dynamics, often referred to as the flow, of X in a 
stated . 

E - Edges, possible transitions from states such as from 
q x to q 2 , edges in Figure 1 are labeled “a” and “b”. 

G tj - Guards, conditions that must be met at an edge to 
transition to another state. 

R t j - Resets, reset maps to define how state variables 
change when transition events from one state to another 
occur, such as from q x to q 2 . 



Reduced ordered models are used here to simplify the 
modeling and control of the biped robot. The sequence 
of motions are divided into action motions, each 
described by a hybrid dynamic system. For example, 
while in a sitting position the biped is modeled as a single 
link, an inverted pendulum. The torso angle is regulated 
to provide correct sitting postures while other actuators 
are left unactuated. When a transition to a standing 
position is desired the controller then attempts to reach 
that desired position needed to propel the biped into a 
desired standing position. The HDS for action motions 
in this paper are described below in Figure 2. 





Figure 2. HDS model of a biped sitting to standing 



Hybrid dynamic subsystems S x and S 2 define sitting and 
standing respectively. Thus 



Q:Q ={q l = sitting } 

X :X I R 2n ,X =[jc p jc 2 ] 



D \D(q l ) = \X / R 2 ,xJ 



3 n n 
~ 2 9 ~2 



> x 2 1 [o,«w] 






</ :/ («„*) = [*]=/, 

E : {(<7i>4i)} 

G :G n = [!s tan d _event] 

R:R n =[x + =x] 



Where x x and x 2 are the angular position and angular 
velocity of the hip, CO Hmax is the maximum angular 

velocity of the hip, and f x are the non-linear equations 
defining the 1-link torso. 



Likewise, the hybrid automaton model S 2 is defined as: 
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Q:Q={q l = sitting, q 2 = standing} 

X : Xj ed 8% ,X 1 = [x,, jc 2 ],n, = 1 

X 2 ^fl 2>h ,X 2 = [x l ,x 2 ,x 3 ,x 4 ],n 2 =2 

D (9 2 ) = | x 2 sO V, 

/:/( ? „X,) = [X 1 ]=/ 1 

/(«2.x 2 )=[xJ=/ 2 
E = {(q l ,q 2 ),{q 2 ,q l )} 

G '■ G n = [x x > 0 m2 \ 

G 21 = < 6 h21 ,x 3 < 0 K2] ] 



R : R n = 



0 

0 

X u 

X,, 



( 2 ) 



■ISI 



Where for states X x ,x x and x 2 are the angular position 
and angular velocity of the hip respectively. For 
states X 2 , X 1 and x 3 are the angular positions of the 

knee and hip respectively, x 2 and x 4 are the angular 

velocity of the knee and hip respectively and CO K max is 
the maximum angular velocity of the knee. The non- 
linear equations in f 2 define the dynamics of the 2-link 

torso and knee. Guard conditions are denoted as 0 Hl2 , 
the hip angle threshold to transition to a standing position, 
and 0 H2l and6^ 21 , the hip and knee angle thresholds to 

be below in order to transition back to a sitting position. 
During a reset, the mapping R defines how states are 
assigned. Where X n and X 12 are the x x and x 2 states 

in the f x dynamics. Similarly, X 23 and X 24 are the x 3 
and x 4 states in the f 2 dynamics. 

3. Biped Robot Model Derivation 

It is well understood that to simulate walking, a biped 
robot model needs five degrees of freedom [3]. Several 
models exist for five link bipeds [l]-[4] and the dynamic 
models while complicated can be described by nonlinear 
differential equations. A general equation describing the 
dynamic motion of a rigid robot takes the form in (3). 

D [q)q +H (q,q)q +G (q ) = a (3) 

D (g) - the inertia matrix of the robot, size of nxn. 

H (q,q) - the friction, coriolis centripetal forces, size 
of nxl. 

G (g) - the gravitational forces acting on each link’ s 
center of gravity, size of nxl. 

u - the actuator torques acting on each link, size of nxl. 

The differential equations can be derived using 
Lagrange’s method for rigid bodies [5]. The equations of 
a system are shown in (4) modified from Zaher [6]. The 

HDS system S x in Figure 2(a) represents a sitting biped 
and thus only the dynamics of the torso are modeled. 



Figure 3 represents the single link torso model. For 
clarity the leg is shown but not annotated. 




Figure 3. Free body diagram of the single link torso 
model 



J. 

2 



m T ljd H +d H 0 H +—m T gl T sin 0 H =u H 



( 4 ) 



Where m T is the mass of the torso, l T is the length of 



the torso, d H is the damping coefficient of the hip joint, 

g is the acceleration due to gravity and 0 H , 0 H , 0 H are 
the angular position, velocity and acceleration of the 
torso respectively. The input U H is the torque applied at 
the hip. A state space form of (4) can be shown with 
states X x = 0 H and x 2 = 0 H . The state space 

representation is used in section 4 to design a 
backstepping controller. Equation (5) also represents the 

dynamics f x in Figure 2(a) in HDS S ] and Figure 2(b) in 



HDS S 2 . 



x l = x 2 

x 2 = -a sin x x - bx 2 + cu 

3 g 3 d H 3 

Where a = , b = — and c = - 

21 T m T l T TYijlj 



( 5 ) 



The HDS system S 2 in Figure 2(b) represents a biped 



beginning to rise from the sitting position to an upright 
standing position. Both right and left legs and feet are 
kept together. Figure 4 shows the double link model of 
the torso and upper leg. As before, the derivation of the 
dynamics is done using Lagrange’s method. The 
equations of this system are shown in equations (6)-(12) 
which are modified from Zaher [6]. 




Figure 4. Free body diagram of the two link torso and upper leg model 
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u K 
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IAJ 
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G„_ 
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Where: 

d n = I UL + I T + (wi T + m uL + ^ulKt cos ) 

d X2 — d 2X — Ij + (jn T + yyi UL ) IjjiJct cos 

22 ^ T 



1 I f vyi _L 9/77 ^ 

__|_r r _ __ 7 iii'T t ±rn UL 
l cUL 0 l UL ’ L cT 0 G 

2 ^ v y 

I UL — ~ m u l JuL ’ I T — ~^ m T^T m UlJl 



( 8 ) 

(9) 



// = (m T + m^ L ) l UL l cT sin 0 H (10) 

= (jn T + m UL ) g/ cr sin + 0 H ) (11) 

= G h + m UL gl cUL sin 0 H + ( m T + ) gl UL sin 0 K (12) 

Where m UL = 1.5kg, m r = 1.5kg 

l UL = 0.3m and l T = 0.2m . Where m i;/ is the mass 

of the upper leg, l UL is the length of the upper leg, d K is 



the damping coefficient of the knee and 0 K , 0 K , 0 K are 
the angular position, velocity and acceleration of the knee 
respectively. The input U K is the torque applied at the 
knee. A state space form of (6) can be shown with 
states x, = 0 K , x 2 = 0 K , x 3 — 0 H , x 4 = 0 H , 

u x — U K and u 2 = U H . The state space representation 
is used in section 4 to design a backstepping controller 
and is also the representation of f 2 in Figure 2(b) in 



HDS S 2 . 



1 2 

1 3 

A. 



x 2 




"0“ 


A 3 + dy 2 x 4 + c 121 x 3 x 4 + c 211 x 4 x 3 + c 221 x 2 + Gj 


+ 


u x 


*4 




0 


d 2l Jc 3 + d 22 Jc 4 + c 112 x^ + G 2 




u 2 



(13) 



Where the Christoffel symbols, terms c ijk in equation(13), 
are determined from the following equations from [5]. 

dd h 



^ijk 



dd kj 



+ - 



dd i} 



dq t dcj j dq k 



(14) 



4. Control In Hybrid States 

Nonlinear backstepping control is used around each 
hybrid state. This method of control involves designing 
virtual stabilizing control functions for each state by 
using Lyapunov functions [18]. Stability and 
performance are considered during design with the 
choice of both the Lyapunov function candidates and the 
control parameters of the stabilizing functions. Stability 
of the overall system is guaranteed since the design 
results in a decomposed sequence of stable Lyapunov 
functions defining the overall system. 



The dynamics of hybrid system S x can be re-written from 
(5) in the following form: 



*1 =*2 



(15) 



x 2 = f(x x ,x 2 ) + gu 

Where f (. 1 ', ,x 2 \ = —a sin x, — bx 2 is a typical 
nonlinear equation describing the 1-link system. The 
first step of the backstepping design is to stabilize the x x 

substate dynamics with a stabilizing function, OC x ( ) . 

We augment the design by adding integral action to help 
reject disturbances and reduce steady state error. This 
design method is similar to the work in [11] . A new 

state variable Z 1 = x ld — X 1 is introduced where x ld is 
the desired hip angle to track. By taking time derivative: 

Zi=X ld ~ X 2 ( 16 ) 

An error variable z 2 is then defined. 

z 2 = a l ()-x 2 (17) 

Substituting (17) into (16) gives the dynamics of the new 
state z x in terms of z 2 : 

z, = x ]d - a, ( ) + z 2 (18) 

Next a control Lyapunov function (elf) is used to 
determine (X x ( ) . The following elf is considered: 

J. 

2 " 2 

The first term takes into account the integral of the 
tracking error, 



v,=^C 2 +-kzj 



(19) 



C = \z,(r)dT . 



( 20 ) 



( 21 ) 



Taking the time derivative of (19) results in: 

Vi =k 1 £z 1 + z x z x 

= z l (k l C + x ld -a l ( ) + z 2 ) 

For the state z 1 to be stable, OC x ( ) must be defined to 
make V x < 0 . Consider the following possible 



stabilizing function where c x and k x are both positive 
control parameters: 

«i( ) = c l Z l +x ld +k l C (22) 

Taking the time derivative of (22) results in: 
a j( ) = c x z x + x Xd + k x z x (23) 

Substituting (22) into (21) results in: 

V t = z, (k x C + x Xd - ( c x z x + x Xd + k x C) + z 2 ) 

= z x (-c x z x +z 2 ) (24) 

= -c x z* + z x z 2 

The state z x is asymptotically stable when the error z 2 
is zero, < 0 . 
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The second step is to stabilize the x 2 sub state dynamics 

using a similar procedure that was used in step 1 . 
Differentiating (17) and from state equations, results in 

the dynamics of the error variable z 2 • 
z 2 = a i ( )-x 2 

, , , v (25) 

= a i\ )- f\ x v x i)-8 u 
The following elf is considered to stabilize the second 
state: 

V 2 =v,+^ (26) 

Taking the time derivative of (26) results in: 

V 2 = Vj + z 2 z 2 

(27) 

= -crf + Z 1 z 2 + z 2 (d 1 ( )-f(x 1 ,x 2 )-gu) 

Since the control input U is part of the error dynamics no 
additional stabilizing function is required. Instead, U is 
defined to stabilize the z 2 state dynamics by 

making V 2 < 0 . This is done by introducing another 



control parameter, c 2 , and by defining the desired error 
dynamics. 

Z 2 = -C 2 Z 2 - Z x (28) 

Then by using (25) and (28) the control input can be 
determined. 



«l( )-f{ x i’X 2 )-gu = -c 2 z 2 -z l (29) 

)~f( X l’ X 2) + C 2Z 2 + Zl) 

u = —(cl sin x l + bx 2 + x u + (k x + 1 — c\ ) z x — + (q + c 2 ) z 2 ) 



Substituting (29) back into (27) results in: 

V 2 = -Cjzf - c 2 zl (30) 

When c 2 is a positive control parameter V 2 < 0 . 

Control parameters C, and c 2 are set to adjust the desired 
behavior of the asymptotically decreasing position and 
velocity errors respectively. The parameter k x sets the 

gain on the added effect resulting from integral action. 

A diagram of the original system of equations (15) is 
shown in Figure 5. The additions to the system diagram 
with backstepping and integral action are shown in 
Figure 6. The control function u for the backstepping 
controller is shown in Figure 7. 




Figure 5. System diagram of the original system S 



Xid %ld 




Figure 6. System diagram showing backstepping with integral action 
control for system S ^ 



xi d 






*1 




Xi 

*2 



Figure 7. Control function for backstepping with integral action for 



system S ^ 

The dynamics for f 2 in the upright state of HDS S 2 can 
be rewritten from (13) in the following form: 



*1 




x 2 




'o' 


x 2 




/*( ) 


4- 


U j 


x 3 




*4 




0 


_*4_ 




M ). 




u 2 



Where f k ( ) and f h ( ) are nonlinear equations 

describing the knee and hip dynamics respectively. The 
first two steps of the backstepping design follow from the 
single link example with some changes to the variable 
subscripts to allow for additional states and control inputs. 

A new state variable z 1 = x ld — X x is introduced where 

X ]d is the desired knee angle to track. By taking time 
derivative: 

Zi=x u -x 2 (32) 

An error variable z 2 is then defined. 

z 2 = a i ( )-x 2 (33) 

Substituting (33) into (32) gives the dynamics of the new 
state z x in terms of z 2 : 



z x =x ld -a x { ) + z 2 (34) 

Next a control Lyapunov function (elf) is used to 
determine OC x ( ) . The following elf is considered: 

v,= l^ + h' (35) 

The first term takes into account the integral of the 
tracking error, 

t 

£i=Jzi {r)dt . (36) 

0 

Taking the time derivative of (35) results in: 
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Vi = k Az x + zA 

, ! \ \ (3?) 
= z 1 {krf x + x xd -a x ( ) + z 2 ) 

For the state z l to be stable, GC X ( ) must be defined to 

make V x < 0 . Consider the following possible 

stabilizing function where C ] and k ] are both positive 
control parameters: 

( ) - c \^\ + X\d + KCi (38) 

Taking the time derivative of (38) results in: 

«1 ( ) = c x z x + X u + k x z x (39) 

Substituting (38) into (37) results in: 

= Zl ( k A + Xu - ( C l^l + Xu + KO + z 2 ) 

= z ] (-c ] z ] +z 2 ) (40) 

= -C,Z 2 + z x z 2 

1 l 1 z 

The state z 1 is asymptotically stable when the error z 2 
is zero, V x < 0 . 

The second step is to stabilize the x 2 sub state dynamics 

using a similar procedure that was used in step 1 . 
Differentiating (33) and from state equations, results in 

the dynamics of the error variable z 2 : 
z 2 =a x ( )-x 2 

= «i( )“/*( )~ u i 

The following elf is considered to stabilize the second 
state: 



(41) 



v 2 =v l +i-g 



1 

2 



(42) 



(43) 



Taking the time derivative of (42) results in: 

y 2 = Vi + z 2 z 2 

= -c x zf + z,z 2 + z 2 (d, ( )-/ t ( )-Mj) 

Since the control input U x is part of the error dynamics 
no additional stabilizing function is required. Instead, U x 
is defined to stabilize the z 2 state dynamics by making 
V 2 < 0 . This is done by introducing another control 
parameter, c 2 and by defining the desired error dynamics. 

z 2 = -C 2 Z 2 - Z x (44) 

Then by using (41) and (44) the control input can be 
determined. 

«l( )-/*( )-«! =-c 2 z 2 -z, 

/X /X (45) 

U l =-/*( ) + «l( ) + C 2Z2+Z 1 
Substituting (45) back into (43) results in: 

V 2 = -c x z x - C 2 z\ (46) 

When c 2 is a positive control parameter V 2 < 0 . 

The following two steps design the backstepping control 
for the states describing the hip dynamics. A new state 



variable z 3 = x 3d — x 3 is introduced where x 3d is the 
desired hip angle to track. By taking time derivative: 



^3 — X 3d X 3 

An error variable z 4 is then defined. 



(47) 



z 4 =a 2 ( )■ 



(48) 

Substituting (48) into (47) gives the dynamics of the new 
state z 3 in terms of z 3 : 

z 3 = x 3d -a 2 ( ) + z 4 (49) 

Next a control Lyapunov function (elf) is used to 
determine a 2 ( ) . The following elf is considered: 

K, r 2 i \ 

2 42 2 

The first term takes into account the integral of the 
tracking error, 



v 3 =v 2 +^!+i^ 



(50) 



i 

£ 2 =\z 3 (T)dT . 



(51) 



(52) 



Taking the time derivative of (50) results in: 

V 3 =V 2 + k 2 C 2 z 3 + z 3 z 3 

= -C x z x - c 2 z 2 + z 3 {kiCi + h d -a 2 ( ) + z 4 ) 

For the state Z 3 to be stable, a 2 ( ) must be defined to 

make V 3 < 0 . Consider the following possible 

stabilizing function where c 3 and k 2 are both positive 
control parameters: 

^2 ( ) — C 3^3 ~l~ X 3d ~l~ k 2 C 2 

Taking the time derivative of (53) results in: 

d 2 ( ) = c 3 z 3 + x 3d +k 2 z 3 (54) 

Substituting (53) into (52) results in: 

V 3 = -erf - erf + Z 3 (krf + x M - (c 3 z 3 + x 3d + krf ) + z 4 ) 

= ~ C rf - C l4 + Z 3 (- C 3^3 + Z 4 ) 



(55) 



CjZi ^2^2 C3Z3 3 ” ^3^4 



The state z 3 is asymptotically stable when the error z 4 
is zero, V 3 < 0 . 

The final step is to stabilize the x 4 sub state dynamics 

using a similar procedure that was used in step 2. 
Differentiating (48) and from state equations, results in 

the dynamics of the error variable z 4 : 
z 4 = d 2 ( ) - x 4 

= « 2 ( )“/*( )~ U 2 

The following elf is considered to stabilize the second 



(56) 



state: 



1 



V =V + — z 2 

V 4 y 3 2 ^ 

Taking the time derivative of (57) results in: 



( 57 ) 
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V 4 = v 3 +^4 

=-c^-c 1 ^-c 3 z 3 +z 3 z A +z 4 (a 1 { )-f h { )-u 2 ) 

Since the control input U 2 is part of the error dynamics 
no additional stabilizing function is required. Instead, 
u 2 is defined to stabilize the z 4 state dynamics by 

making V 4 < 0 . This is done by introducing another 

control parameter, c 4 and by defining the desired error 
dynamics. 

Z 4 = -C 4 Z 4 - Z 3 (59) 

Then by using (56) and (59) the control input can be 
determined. 

« 2 ( )"/*( )- U 2=-C^4~Z Z 

U 2 = ~fh ( ) + ^2 ( ) + C 4^4 + Z 3 

Substituting (60) back into (58) results in: 

V 4 = —C X Z X — C 2 Z 2 — C3Z3 — c 4 z 4 
When c 4 is a positive control parameter V 4 < 0 . 

Control parameters C, and c 2 are set to adjust the desired 

behavior of the asymptotically decreasing position and 
velocity errors respectively of the knee. Likewise the 

control parameters c 3 and c 4 are for the hip. The 

parameters k x and k 2 set the gain on the added effect 

resulting from integral action for the knee and hip 
respectively. 



(60) 

(61) 



5. Simulations 

Sitting to standing posture simulation is performed with 
the sitting portion over a fixed period of time followed by 
a standing portion with another fixed period of time. 
Each fixed portion of time represents the time in HDS 

state S ] and S 2 with a standing event triggered between 

the fixed time intervals. During each HDS state a 
disturbance is introduced to examine robustness. This 
disturbance represents a force applied on the biped torso 

during state S x and applied to the upper leg during 

state S 2 . Robustness to system parameter variations and 
the benefits of integral action are also demonstrated. 



The reference positions are pre-filtered through a model 
reference to provide the positions and derivatives 
required by the backstepping controller. A second order 
reference model is used to provide realistic trajectories to 
track without overshoots [5]. The use of a second-order 
model also makes it possible to obtain the velocity and 
acceleration information from the reference without the 
need to specify these trajectories. The following model 
was used from Khalil 



[18]. 



Ti" 




y 2 


_f 2 _ 







(62) 



Parameter of the model are set to CO n =10 and £ = 1 to 

give a critically damped response with approximately a 
0.5s settling time. 

Initial positions for S x are 0 H (0) = 0 , 

0 K ( 0) = — — which simulates a sitting position with 
the torso completely vertical. The desired position for 

5, is c/ w = — which is a slightly leaned forward torso. 

des 6 

Once the transition from state S x to S 2 occurs the desired 

/i 71 

torso position changes to u H = — which is a deeply 

leaned forward position to get the torso center of gravity 
above the position of the foot. This position is in 

preparation to stand upright. Once S 2 has transitioned to 
the upright sub-state, f 2 , the desired positions change to 
the final upright standing posture of 0 H = 0 , 



Actuator Torque -Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 8. 1-link backstepping response for HDS S x 

Figure 8. shows the response during the state S x . The 
actuator torque, desired and actual positions, angular 
velocity and a phase plane plot of 0 H vs. 0 H is also 

shown. From the angular tracking response plot is seen 
that the model reference trajectory is followed closely. 
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At t = 2s a 3Nm disturbance torque is applied for 
200ms. The resulting correction in tracking can be seen 
as the loop in the phase plot. 



Actuator Torque - Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 9. 1-link backstepping response for HDS S 2 f\ 

Figure 9. shows the response during the state S 2 when 

the dynamics j\ are present. This is the hip only portion 

of motion where the desired position is a deep lean 
forward. From the angular tracking response plot it is 
seen that the model reference trajectory is followed by 
the hip. 

Figure 12. shows the total system response with four 
second fixed periods of time in states 5 1 and S 2 . A 

standing event is simulated at t = 4s . Disturbance 
pulses of 3Nm are present for 200ms at t = 2s on the 
torso and at t = 6s on the upper leg. The reference and 
actual position of the hip and knee along with actuator 
torque is shown. Notice that during state S l the knee is 

assumed fixed in position since only the torso is moving 
on the seated biped. No actuator torque is applied to the 

knee until in state S 2 when in the upright sub-state, q 2 of 
Figure 2(b). 



Actuator Torque - Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 10. 2-link backstepping response of the hip for HDS S 2 f 2 



Actuator Torque - Knee 




Angular Tracking Response - Knee 




Angular Velocity - Knee 




Figure 11. 2-link backstepping response of the knee for HDS S 2 f 2 
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Angular Tracking Response - Knee 




Angular Tracking Response - Hip 




Actuator Torque - Knee and Hip 




Figure 12. Total system response over HDS S ^ to S 2 using 
backstepping control 



This simulation uses the following controller parameters 

Cj,c 3 ,c 3 ,c 4 = 3 ,k v k 2 =4 ,a = 29.43, b = 3,c = 6 

Figure 13. shows the same system simulated with PD 
control instead of using the backstepping control from 
section 4. The backstepping control is much smoother in 
following the references trajectories when compared to 
the PD control. Additionally, there is less output 
chattering and lower output magnitudes required when 
using the backstepping control. 

Simulink and Simmechanics diagrams of the system are 
shown in figure 14 and figure 15 respectively. 
Simmechanics was established to be a representative 
method to simulate robotic systems by Dung and Kang 
[14]. The authors have also shown the use of 
Simmechanics in simulations of hybrid dynamic systems 
in their previous work [19]. 





Actuator Torque - Knee and Hip 




Time(sec.) 



Figure 13. Total system response over HDS to using PD control 




Enijifronnieitf 



Figure 14. Simulink diagram of the simulated system 




Figure 15. Simmechanics model of biped 
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5.1 Robustness 

Simulations to demonstrate robustness are divided into 
two parts. First, the robustness to parameter variations of 
the system is examined. Second, the robustness to 
disturbance is examined. Parameters for the mass, length 
and joint damping are changed +/- 30% to see the effects 

on the response. Phase plots during state S ] are shown in 

the following figures, Figure 16-18, to illustrate the 
response differences as parameters are adjusted. 



Atfuator Torque -Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 16. Robustness to variations in system mass 



The controller is robust to mass variation in the system. 
A system with less mass responds quicker and therefore 
is more difficult to control with the same actuator. 
However, the response does show that even for the -30% 
mass variation the system is asymptotically stable. 

Robustness as a result of length variation is very similar 
to mass variations. A shorter link responds quicker and 
is more difficult to control with a given actuator. Again, 
the response for the reduced length variation is still 
asymptotically stable. 



Actuator Torque -Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 17. Robustness to variations in system length 



Actuator Torque - Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 18. Robustness to variations in system damping 
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Robustness to damping variation is the best of the three 
parameters. The effect of damping does not contribute as 
much to the system response and therefore the controller 
is able to compensate changes well. 

5.2 Disturbance Rejection 

Disturbance during the SI and S2 states is simulated by a 
step function of applied toque. The simulations are 
repeated with and without the integral action of the 
backstepping controller to show the performance 
difference. The disturbance signal is shown in Figure 19. 
The amplitude is 3Nm and the duration is 200ms. 
Trajectory tracking performance is also examined to see 
the benefits of the integral action. 

3Nm 



t = 200 ms t 

Figure 19. Disturbance signal 



Actuator Torque - Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




The integral action improves the response time to reject 
the disturbance effect with the consequence of adding 
overshoot. This is seen in Figure 20(b). 

Figure 21 examines the trajectory tracking performance 
of the integral action. A trapezoidal trajectory is used 
where there is a constant slope from t = 0s to t = 2s . 
Then at t = 2s the reference is at a fixed value of 
71 

— radians. Again as in Figure 20(b) the response with 

6 

integral action is faster with some overshoot. The error 
and sum squared error is examined further in Figure 22. 
The error between the responses with and without 
integral action appears to show that without integral 
action is better. However, looking at the sum square 
error the integral action response does perform better 
during the ramp portion of the trajectory tracking and has 
an overall smaller sum square error. 



Actuator Torque -Hip 




Angular Tracking Response - Hip 




Angular Velocity - Hip 




Figure 21. Response comparing trajectory tracking with and without 
integral action 



Figure 20. Response comparing disturbance rejection with and without 
integral action 
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ErrorTracking Response - Hip 




Sum Squared ErrorTracking Response - Hip 




Figure 22. Error and sum squared error of trajectory tracking with and 
without integral action 

6. Conclusion 

Sequenced hybrid dynamic systems have been shown to 
be useful in the control of simulated biped robot motion 
from a sitting to standing pose. Backstepping control 
with integral action was applied to control the biped in 
the individual hybrid states. Robustness of the 
backstepping control with integral action was shown to 
be robust to variations in parameters in mass, length and 
damping with an asymptotically stable response. The 
benefit of integral action to the backstepping control was 
shown to improve the response time and results in an 
overall better performance during trajectory tracking with 
regards to the sum square error of the response. Overall 
the control method presented when compared to a 
classical PD control is better from trajectory tracking and 
smoother output actuator torque. Extensions to this work 
with additional sequenced poses and motions are planned. 
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